Ricatti方程式を解く

[縁の下の力持ち]

初めまして。規約を読んでどう書いたらいいのか分からなかったのでコピペして下さい。というものをコピーして質問したいと思います。

[1] 質問文
　[1.1] 自分が今行いたい事は何か
　Ricatti方程式を数値的求解法を用いて安定化ゲインFを設計することです。
　[1.2] どのように取り組んだか(プログラムコードがある場合記載)
　C言語については全く分からない状態です。
　[1.3] どのようなエラーやトラブルで困っているか(エラーメッセージが解る場合は記載)
　他学科履修で授業を受講したため、C言語について何も分からない状態で課題が出て困っています。
　[1.4] 今何がわからないのか、知りたいのか
　C言語については全く分からない状態なのでお手上げ状態です。授業の際にレポートについてヒントを先生が　言っていて、for文？とかを使うそうです。曖昧で申し訳ありません・・・

[2] 環境 　
　[2.1]OS ： Windows XP
　[2.2]コンパイラ名 ： Borand C++を使うそうです。
[3] その他
　・どの程度Ｃ言語を理解しているか
　　情報系の学科出身ではないため、C言語については全く分からない状態です。何かヒントなどを頂ければこの機会に学ぼうと思います。

丸投げは禁止とあったのですが、C言語について何も分からない状態のため投稿させて頂きました。もしこの掲示板に不適切であれば管理人様削除して頂いて結構です。

もちろん答えを示して頂ければ幸いですが、もし宜しければC言語について何か教えて頂ければ嬉しいです。宜しくお願いいたします。

[縁の下の力持ち]

追記です。

[TOMONORI]

すみません。リッカチはともかく、安定化ゲインてなんですか？

ググってみたらオクターブでの解はそれっぽいのがありました。
参考になるかは不明ですが・・・
ttp://anabuki.ec.u-tokai.ac.jp/class/class/2006/system/system_1.pdf
１９ページ（２６枚目）辺り。

[縁の下の力持ち]

レスありがとうございます。
安定化ゲインFはおそらく追記の画像で示した
J=～の評価関数を最小にする値だと思います。

[たいちう]

当面の問題ははっきりしていて「C言語について何も分からない」ことですね。
まずそこを何とかしましょう。

入門書やサイトがいくらでもありますので、それを勉強してください。
C言語全般を勉強することが理想ですが、
とりあえずは二次方程式を解けるくらいまでで十分かと。
二次方程式のC言語のプログラムなど、Web上で簡単に見つけられると思いますが、
見つけろといっているのではないですよ。自分で理解して解法をプログラムしてください。

初心者用のHPを1つ紹介します。自分にあったものを見つけてください。
http://www9.plala.or.jp/sgwr-t/index.html

[縁の下の力持ち]

>>たいちうさん

参考サイトのURLありがとうございます。
この機会にC言語を学ぼうと思います。
ただ、期限までにプログラムが作成できるまで身につくか心配ですが・・・

[non]

コンピュータは、ご存じのように、非常に高速に四則演算や論理演算を行う機械です。
人間がやることを、代わりに高速におこなってくれるものです。
従って、アルゴリズムは自分で考えなくてはいけません。
実際の課題をあなたは、紙と鉛筆を使ってどう解きますか？
それを示していただければ、それに沿ったプログラムを作るお手伝いをします。

まして、今回のこの課題は、C言語を履修していない学生がいるのにも関わらず出している
のですから、プログラムはなくてもアルゴリズムを提出すればいいのではないのでしょうか？

[縁の下の力持ち]

>>nonさん

コメントありがとうございます。
早速手で解いてみようと思います。そうですよね、手で解けないことを機械の力に頼る感じですよね。
ただ、アルゴリズムの提出だけでは不十分です。
この授業は他学科履修の授業だから受講している学生の中でC言語を知らないのは僕だけでしょうし、この授業を取らないといけなくなったのも自分の責任ですから(汗)

[縁の下の力持ち]

こんな感じで解いてみました。
合ってるか分かりませんが(汗)

これをプログラムで組むとどのようになるのでしょうか？

[non]

制御理論についてはさっぱりわかりませんので少し質問させてください。

１　step１のｆの初期値はどうすればいいのですか？

２　step２の計算途中で、Ｑがなくなっているのはあってますか？

３　step３の収束条件の||P(i)-P(i-1)||ってのは具体的にどういう計算のことですか？
　　特に初期値、i=0のときの式を示してもらえますか。

４　絶対に収束はしますか？

[縁の下の力持ち]

>>nonさん

少し間違っていたので修正した添付ファイルの載せておきます・・・

質問の答えに関してもそこに記してありますが

１　fの初期値は適当に決めるそうです。
２　Qがなくなっていたのは単純に自分の計算ミスです。すいません・・・
３　おそらくP(i)-P(i-1)の差の行列式の値→→→|P(i)-P(i-1)|のことだと思います。
４　収束は絶対にすると思います。

[non]

だいぶ意味がわかってきました。
step２の最後の式は、
ｐ＝の式にして欲しいですね。
たぶん逆行列を左からかけるのでしょうけど。。。
逆行列を変換した式になおしてもらえますか？
なぜっていうと、逆行列をプログラムで求めることはできますが、（掃き出し法とか）
今回は２行２列と決まっているので、計算した式で求めた方が簡単ですから。

[non]

>今回は２行２列と決まっているので
４行４列でしたね。

[縁の下の力持ち]

>>nonさん

すいません。所用で少し離れていました。
4行4列の逆行列を求めて、ｐ＝の式にしたいと思います。

[non]

逆行列を掃き出し法で作っていたのですが、どうもうまくいきません。
誤差が生じているのでしょうかね。
別の方法で逆行列を求める必要がありますが、私に時間がありません。
直接値を入れて逆行列を求めるか、誰か時間のある人、手直ししてください。

[縁の下の力持ち]

>>nonさん

プログラムを作って頂きありがとうございます。
しかし、Borlandで実行してみたところ行列が表示されるだけでゲインfの値が表示されないのですが・・・

もしかしたら誤差はあるのかもしれません。明日聞きに行ってみます。

[non]

上にも書いたように、逆行列が求まっていないのです。
直接、逆行列の計算式を作っていれてもらえますか。

[縁の下の力持ち]

>>nonさん

了解しました。
逆行列を求めてみます。

[non]

同じ内容の行があるので、逆行列が求まらないのでは？
ただ、未知数が3個で、式が３つなので答えは求まるはず。
P11，P12,P22を直接求めて式を入れるのかなぁ？
お友達はどうしてます？

[縁の下の力持ち]

>>nonさん

先ほど友達に確認すると3行3列で計算していました。

[縁の下の力持ち]

3行3列でやったらできました。
一応プログラム載せておきます。

#include <stdio.h>

double A[2][2] = {{-2, 0}, {2, 1}};   //行列A
double B[2] = {-1, 1};                //行列B
double Q[2][2] = {{200, 0}, {0, 200}};  //重み行列Q
double R = 10;                          //重み行列R
double f[2] = {-4, -6};               //設計ゲインF
double Af[2][2];                      //行列AF(i)
double P[2][2];                       //リアプノフ方程式の解P(i)
double oldP[2][2] = {0};              //収束条件計算用P(i - 1)

/* 行列表示関数 */
void disp_array(char *str, double *a, int n, int m)
{
    int i;  //行
    int j;  //列
    
    puts(str);
    for(i = 0; i < n; i++){
        for(j = 0; j < m; j++){
            printf("%f ",*(a + i * m + j));
        }
        printf("\n");
    }
}

/* 行列AF(i)を求める関数 */
void setAf(void)
{
    int i;  //行
    int j;  //列
    
    for(i = 0; i < 2; i++){
        for(j = 0; j < 2; j++){
            Af[j] = A[j] + B * f[j];
        }
    }
}

/* 逆行列を求める関数 */
void inv(double a[/url][3])
{    
    double inv_a[3][3];  //ここに逆行列が入る
    double buf;          //一時的なデータを蓄える
    int i;
    int j;
    int k;
    
    /* 単位行列を作る */
    for(i = 0; i < 3; i++){
        for(j = 0; j < 3; j++){
            if(i == j){
                inv_a[j] = 1.0;
            }else{
                inv_a[j] = 0.0;
            }
        }
    }
    
    disp_array("inv_a:", (double *)inv_a, 3, 3);
    
    /* 掃き出し法 */
    for(i = 0; i < 3; i++){
        buf = 1 / a;
        for(j = 0; j < 3; j++){
            a[j] *= buf;
            inv_a[j] *= buf;
        }
        for(j = 0; j < 3; j++){
            if(i != j){
                buf = a[j];
                for(k = 0; k < 3; k++){
                    a[j][k] -= a[i][k] * buf;
                    inv_a[j][k] -= inv_a[i][k] * buf;
                }
            }
        }
    }
        
    /* 求めた逆行列を代入 */
    for(i = 0; i < 3; i++)
        for(j = 0; j < 3; j++)
            a[i][j] = inv_a[i][j];
}

/* リアプノフ方程式の解P(i)を求める関数 */
void makeP(void)
{
    double a[3][3] = {0};
    double c[3];
    double p[3] = {0};
    int i;
    int j;
    
    a[0][0] = 2 * Af[0][0];
    a[0][1] = 2 * Af[1][0];
    a[1][0] = Af[0][1];
    a[1][1] = Af[0][0] + Af[1][1];
    a[1][2] = Af[1][0];
    a[2][1] = 2 * Af[0][1];
    a[2][2] = 2 * Af[1][1];
    
    disp_array("a:", (double *)a, 3, 3);     //debug
    
    /* 逆行列を計算 */
    inv(a);
    
    disp_array("inv_a:", (double *)a, 3, 3); //debug
    
    c[0] = R * f[0] * f[0] - Q[0][0];
    c[1] = R * f[0] * f[1];
    c[2] = R * f[1] * f[1] - Q[1][1];
    
    /* 行列P(i)を計算する */
    for(i = 0; i < 3; i++){
        for(j = 0; j < 3; j++){
            p[i] += (a[i][j] * c[j]);
        }
    }
    
    P[0][0] = p[0];
    P[0][1] = p[1];
    P[1][0] = p[1];
    P[1][1] = p[2];
}

/* 収束条件計算関数 */
double e_range(void)
{
    int i;           //行
    int j;           //列
    double e[2][2];  //P(i)とP(i - 1)の差を格納する行列
    double ret;      //行列式
    
    // 差を求める
    for(i = 0; i < 2; i++){
        for(j = 0; j < 2; j++){
            e[i][j] = P[i][j] - oldP[i][j];
        }
    }
    
    //行列式
    ret = e[0][0] * e[1][1] - e[0][1] * e[1][0];
    
    //oldP(P(i - 1))へコピー
    for(i = 0; i < 2; i++){
        for(j = 0; j < 2; j++){
            oldP[i][j] = P[i][j];
        }
    }
    
    if(ret < 0)
        return -1 * ret;
    return ret;
}

/* ゲインF(i + 1)を求める関数 */
void makenewf(void)
{
    f[0] = -1 / R * (B[0] * P[0][0] + B[1] * P[0][1]);
    f[1] = -1 / R * (B[0] * P[0][1] + B[1] * P[1][1]);
}

/* 行列AF(i + 1)を求める関数 */
void makenewAf(void)
{
    Af[0][0] = A[0][0] + B[0] * f[0];
    Af[0][1] = A[0][1] + B[0] * f[1];
    Af[1][0] = A[1][0] + B[1] * f[0];
    Af[0][0] = A[1][1] + B[1] * f[1];
}

/* ゲインFを逐次表示する関数 */
void dispf(int i,double e)
{
    printf("%d回目f1 = %f f2 = %f e = %f\n", i, f[0], f[1], e);
}

/* メイン関数 */
int main(void)
{
    int i = 0;
    int ch;
    double e;
    
    setAf();  //AF0を求める。
    disp_array("Af:", (double *)Af, 2, 2);
    
    /* 収束条件を満たすまで繰り返す */
    while(1){
        makeP();
        e = e_range();
        dispf(i, e);
        if(e < 1.0e-6)
            break;
        i++;
        makenewf();
        makenewAf();
    }
    return 0;
}

[non]

おめでとうございます。
コメントも入れられて、若干変更されているところを見ると、プログラムは理解されたようで・・・
それでは、私は、お役ご免と言うことで。。。。

なお、逆行列を求めている部分（掃き出し法）で、行の入れ替えを行っておりません
入れ替えを行うのなら例えば、

void swap( double *wk1, double *wk2)  /* データの入れ替え関数  */
{
    double     w;
    w= *wk1; *wk1= *wk2; *wk2=w;
}
void inv(double a[/url][3])
{    
    double inv_a[3][3]; //ここに逆行列が入る
    double buf,ak; //一時的なデータを蓄える
    int i,j,k,ik; 
    //単位行列を作る
    for(i=0;i<3;i++)
        for(j=0;j<3;j++)
            inv_a[j]=(i==j)?1.0:0.0;
    //掃き出し法
    for(i=0;i<3;i++){
        ak=fabs(a); //対角要素の絶対値
        ik=i;
        for(j=i;j<3;j++){  //絶対値が最大の要素の行を探す
            if(ak<fabs(a[j])){
                ik=j;
                ak=fabs(a[j]);
            }
        }
        if(ik!=i){ //行の入れ替え
            for(j=0;j<3;j++){
                swap(&a[j],&a[ik][j]);
                swap(&inv_a[j],&inv_a[ik][j]);
            }
        }
        //ここからは前と同じ
        buf=1/a;
        for(j=0;j<3;j++){
            a[j]*=buf;
            inv_a[i][j]*=buf;
        }
        //以下略

のようにしてください。ただし、テストしてません。

[縁の下の力持ち]

>>nonさん

ど素人の僕に丁寧に教えて下さりありがとうございます。しかもプログラムまで書いて頂いて・・・何とお礼を言っていいかわかりません。

本当にありがとうございました。