NaNが出る原因が何か分からない

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら  ※コードを貼り付ける場合は [code][/code] で囲って下さい。詳しくはこちら
yukikirin

NaNが出る原因が何か分からない

#1

投稿記事 by yukikirin » 3年前

現在、とあるロボットのシミュレーションを行いたく、コードを書いているのですが、出力が無限大に発散してしまい、NaNが出てしまいます。

細かいロボットの理論が合っているか聞きたいわけではなく、エラーが起きてそうな記述があれば指摘していただきたいと思っています。

例えば、この記述は0で割ってしまう可能性がある記述であるとか、些細な書き間違いなどです。

些細なことでも構いませんので、宜しくお願いします。

コード:

#include <stdio.h>
#include <math.h>

double dob (double,double,double,double,double,double,double,double,double,double);

//PDゲイン
double kp = 60; //位置ゲイン 
double kv = 12; //速度ゲイン 
     
//運動モデル
const double dt = 0.01;    //サンプリングタイム
const double kt = 33.0;     //推力定数
const double omega = 500;   //カットオフ周波数 100

double tor_fb[3] = {0,0,0}; //DOB feedback
double tor4_old[3] = {0,0,0};
double tor5_old[3] = {0,0,0};


//エンドパートパラメータ
double L_AC = 0.07;
double L_BE = 0.1;
double L_HF = 0.1;
double L_HG = 0.07;
double L_OH = 0.13;

double MC = 0.05;
double ME = 0.07;
double MF = 0.05;
double MG = 0.07;
double MH = 0.07;





int main(void) {
    
    file = fopen("/Users/tamuragou/Desktop/sample.txt","w");
    
    double t;                             //時間
    double g = 9.8;                       //重力加速度
    double pie = 3.1416;                  //π
    
    double xw_cmd, yw_cmd, zw_cmd;
    double dxw_cmd, dyw_cmd, dzw_cmd;;
    
    double ddxw_ref, ddyw_ref, ddzw_ref;
    
    double xw_res = 0, yw_res = -0.00707, zw_res = 0.137;
    double dxw_res = 0, dyw_res = 0, dzw_res =0;
    
    double j[3][3];
    double j_inv[3][3];
    
    double d = 0.01, the = pie/4, phi = 0;                      //モータ角度(初期値)
    double delta;
    
    double ddd_ref, ddthe_ref, ddphi_ref;
    
    double A,B;
    
    double m[3];
    double h_tor[3];
    double g_tor[3];
    
    double    d_res = 0,     the_res = 0,    phi_res = 0;
    double   dd_res = 0,    dthe_res = 0,   dphi_res = 0;
    double  ddd_res = 0,   ddthe_res = 0,  ddphi_res = 0;
    
    double    d_res_old = 0,     the_res_old = 0,    phi_res_old = 0;
    double   dd_res_old = 0,    dthe_res_old = 0,   dphi_res_old = 0;
    double  ddd_res_old = 0,   ddthe_res_old = 0,  ddphi_res_old = 0;
    
    double tor_1[3];
    double tor_2[3];
    double tor_3[3];
    
    double ia[3];
    
    double t_dis[3];
    double t_dis_hat[3];
    
    double ia_cmp[3] = {0,0,0};
    
    
    
    
    
    for (t = 0;t <= 2;t += dt){
        
        
        
        //Target position

            xw_cmd = 0.0001*t;
            yw_cmd = -0.01*sin(pie/4 + 0.01*t);
            zw_cmd = 0.13 + 0.01*cos(pie/4 + 0.01*t);

            dxw_cmd =  0.0001;
            dyw_cmd = -0.001*cos(0.785398 + 0.01*t);
            dzw_cmd = -0.001*sin(0.785398 + 0.01*t);

        
        
        
        //controller(PD制御)
        ddxw_ref = kp*(xw_cmd - xw_res)+kv*(dxw_cmd - dxw_res);
        ddyw_ref = kp*(yw_cmd - yw_res)+kv*(dyw_cmd - dyw_res);
        ddzw_ref = kp*(zw_cmd - zw_res)+kv*(dzw_cmd - dzw_res);
        
        
        //Forward Kinematics
        //ヤコビアンJ
        j[0][0] = sin(the)*sin(phi);
        j[0][1] = d*cos(the)*sin(phi);
        j[0][2] = d*sin(the)*cos(phi);
        
        j[1][0] = -sin(the)*cos(phi);
        j[1][1] = -d*cos(the)*cos(phi);
        j[1][2] =  d*sin(the)*sin(phi);
        
        j[2][0] = cos(the);
        j[2][1] = -d*sin(the);
        j[2][2] = 0;
        
        
        //Inverse Kinematics
        delta = j[0][1]*j[1][2]*j[2][0] + j[0][2]*j[1][0]*j[2][1] - j[0][2]*j[1][1]*j[2][0] - j[0][0]*j[1][2]*j[2][1];
        
        j_inv[0][0] = ( -j[1][2]*j[2][1] ) / delta;
        j_inv[0][1] = ( j[0][2]*j[2][1] )  / delta;
        j_inv[0][2] = ( j[0][1]*j[1][2] - j[0][2]*j[1][1] ) /delta;
        
        j_inv[1][0] =  ( j[1][2]*j[2][0] )  / delta;
        j_inv[1][1] =  ( -j[0][2]*j[2][0] ) / delta;
        j_inv[1][2] = -( j[0][0]*j[1][2] - j[0][2]*j[1][0] ) /delta;
        
        j_inv[2][0] =  ( j[1][0]*j[2][1] - j[1][1]*j[2][0] ) /delta;
        j_inv[2][1] = -( j[0][0]*j[2][1] - j[0][1]*j[2][0] ) /delta;
        j_inv[2][2] =  ( j[0][0]*j[1][1] - j[0][1]*j[1][0] ) /delta;
        
        
        //θ” = J^-1 * x"
        ddd_ref   = j_inv[0][0]*ddxw_ref + j_inv[0][1]*ddyw_ref + j_inv[0][2]*ddzw_ref;
        ddthe_ref = j_inv[1][0]*ddxw_ref + j_inv[1][1]*ddyw_ref + j_inv[1][2]*ddzw_ref;
        ddphi_ref = j_inv[2][0]*ddxw_ref + j_inv[2][1]*ddyw_ref + j_inv[2][2]*ddzw_ref;
        
        
        A = MC*L_AC*L_AC + ME*L_BE*L_BE + MF*L_HF*L_HF + MG*L_HG*L_HG;
        B = MC*L_AC + ME*L_BE + MF*L_HF + MG*L_HG;
        
        
        //Inverse dynamics T = Mθ” + H + G
        //M
        m[0] = MH;
        m[1] = A + MH*d*d;
        m[2] = ( A + MH*d*d )*sin(the)*sin(the);
        
        //H
        h_tor[0] = -MH*d*( ( dthe_res*dthe_res ) + ( sin(the)*sin(the) )*( dphi_res*dphi_res ) );
        h_tor[1] = -( A + MH*d*d )*sin(the)*cos(the)*( dphi_res*dphi_res );
        h_tor[2] = 0;
        
        //G
        g_tor[0] = -MH*g*sin(the)*cos(phi);
        g_tor[1] = ( B - MH*d )*g*cos(the)*cos(phi);
        g_tor[2] = -(B - MH*d )*g*sin(the)*sin(phi);
        
        
        
        
        //T1 = M*θ”ref
        tor_1[0] = m[0]*ddd_ref;
        tor_1[1] = m[1]*ddthe_ref;
        tor_1[2] = m[2]*ddphi_ref;
        
        //Ia = T1/Kt
        ia[0] = tor_1[0]/kt;
        ia[1] = tor_1[1]/kt;
        ia[2] = tor_1[2]/kt;
        
        
        //外乱補償電流を加える
        ia[0] += ia_cmp[0];
        ia[1] += ia_cmp[1];
        ia[2] += ia_cmp[2];
        
        
        //T2 = Kt*Ia
        tor_2[0] = kt*ia[0];
        tor_2[1] = kt*ia[1];
        tor_2[2] = kt*ia[2];
        
        //外乱 Tdis = H + G
        t_dis[0] = h_tor[0] + g_tor[0];
        t_dis[1] = h_tor[1] + g_tor[1];
        t_dis[2] = h_tor[2] + g_tor[2];
        
        
        //T3 = T2 + Tdis
        tor_3[0] = tor_2[0] + t_dis[0];
        tor_3[1] = tor_2[1] + t_dis[1];
        tor_3[2] = tor_2[2] + t_dis[2];
        
        
        //θ"res = T3/M
        ddd_res   = tor_3[0]/m[0];
        ddthe_res = tor_3[1]/m[1];
        ddphi_res = tor_3[2]/m[2];
        
        //θres = 1/s * θ'res
        d_res   = d_res_old   + (dd_res   + dd_res_old)*dt/2.0;
        the_res = the_res_old + (dthe_res + dthe_res_old)*dt/2.0;
        phi_res = phi_res_old + (dphi_res + dphi_res_old)*dt/2.0;
        
        
        //θ'res = 1/s * θ"res
        dd_res   = dd_res_old   + (ddd_res   + ddd_res_old)*dt/2.0;
        dthe_res = dthe_res_old + (ddthe_res + ddthe_res_old)*dt/2.0;
        dphi_res = dphi_res_old + (ddphi_res + ddphi_res_old)*dt/2.0;
       
        
        //前回値としてθ"resを記憶
        ddd_res_old   = ddd_res;
        ddthe_res_old = ddthe_res;
        ddphi_res_old = ddphi_res;
        
        //前回値としてθ'resを記憶
        dd_res_old   = dd_res;
        dthe_res_old = dthe_res;
        dphi_res_old = dphi_res;
        
        //前回値としてθresを記憶
        d_res_old   = d_res;
        the_res_old = the_res;
        phi_res_old = phi_res;
        
        
        //Tdis^ (DOB)
        t_dis_hat[0] = dob(ia[0], ia[1], ia[2], dd_res, dthe_res, dphi_res, m[0], m[1], m[2], 0);
        t_dis_hat[1] = dob(ia[0], ia[1], ia[2], dd_res, dthe_res, dphi_res, m[0], m[1], m[2], 1);
        t_dis_hat[2] = dob(ia[0], ia[1], ia[2], dd_res, dthe_res, dphi_res, m[0], m[1], m[2], 2);
        
        
        //外乱補償電流Iacmp
        ia_cmp[0] = (1/kt) * t_dis_hat[0];
        ia_cmp[1] = (1/kt) * t_dis_hat[1];
        ia_cmp[2] = (1/kt) * t_dis_hat[2];
        
        
        //x'res = J * θ'res
        dxw_res = j[0][0]*dd_res + j[0][1]*dthe_res + j[0][2]*dphi_res;
        dyw_res = j[1][0]*dd_res + j[1][1]*dthe_res + j[1][2]*dphi_res;
        dzw_res = j[2][0]*dd_res + j[2][1]*dthe_res + j[2][2]*dphi_res;
        
        
        //xresをθresから求める
        xw_res =  d_res*sin(the_res)*sin(phi_res);
        yw_res = -d_res*sin(the_res)*cos(phi_res);
        zw_res = L_OH + d_res*cos(the_res);
        
        
        d = d_res;
        the = the_res;
        phi = phi_res;

        
        
//        printf("%lf\t%lf \n",t,zw_res);
//        fprintf(file,"%lf\t%lf  \n",t,zw_res);
    }
    
}


//外乱オブザーバー
double dob (double ia_dob0, double ia_dob1, double ia_dob2, double d_dob, double the_dob, double phi_dob, double m0, double m1, double m2, double count){
    
    double tor_d1[3];
    double tor_d2[3];
    double tor_d3[3];
    double tor_d4[3];
    double tor_d5[3];
    double tdis_dob[3];
    
    
    
    //T1 = kt*Ia
    tor_d1[0] = kt*ia_dob0;
    tor_d1[1] = kt*ia_dob1;
    tor_d1[2] = kt*ia_dob2;
    
    
    //T2 = ωc*(MN * θ'res)
    tor_d2[0] = omega*( m0*d_dob );
    tor_d2[1] = omega*( m1*the_dob );
    tor_d2[2] = omega*( m2*phi_dob );
    
    
    //T3 = T1 + T2
    tor_d3[0] = tor_d1[0] + tor_d2[0];
    tor_d3[1] = tor_d1[1] + tor_d2[1];
    tor_d3[2] = tor_d1[2] + tor_d2[2];
    
    
    //T4 = ωc*T3 - T_feedback
    tor_d4[0] = ( omega*tor_d3[0] ) - tor_fb[0];
    tor_d4[1] = ( omega*tor_d3[1] ) - tor_fb[1];
    tor_d4[2] = ( omega*tor_d3[2] ) - tor_fb[2];
    
    
    //LPF///////////////////////////////////////
    //T5 = 1/s * T4
    tor_d5[0] = tor5_old[0] + ( tor_d4[0] + tor4_old[0] )*dt/2.0;
    tor_d5[1] = tor5_old[1] + ( tor_d4[1] + tor4_old[1] )*dt/2.0;
    tor_d5[2] = tor5_old[2] + ( tor_d4[2] + tor4_old[2] )*dt/2.0;
    
    
    //T_feedback = ωc * T5 (フィードバック)
    tor_fb[0] = omega*tor_d5[0];
    tor_fb[1] = omega*tor_d5[1];
    tor_fb[2] = omega*tor_d5[2];
    
    
    //前回値としてT4を記憶
    tor4_old[0] = tor_d4[0];
    tor4_old[1] = tor_d4[1];
    tor4_old[2] = tor_d4[2];
    
    
    //前回値としてT5を記憶
    tor5_old[0] = tor_d5[0];
    tor5_old[1] = tor_d5[1];
    tor5_old[2] = tor_d5[2];
    
    //.//////////////////////////////////////////
    
    
    //Tdis = T5 - T2
    tdis_dob[0] =  tor_d5[0] - tor_d2[0];
    tdis_dob[1] =  tor_d5[1] - tor_d2[1];
    tdis_dob[2] =  tor_d5[2] - tor_d2[2];
    
    
    //出力
       if(count == 0){
           return tdis_dob[0];
           
       }else if (count == 1){
           return tdis_dob[1];
           
       }else if (count == 2){                     //default
           return tdis_dob[2];
           
       }else{                 //default
           return tdis_dob[2];
       }
    
}
    


アバター
あたっしゅ
記事: 664
登録日時: 13年前
住所: 東京23区
連絡を取る:

Re: NaNが出る原因が何か分からない

#2

投稿記事 by あたっしゅ » 3年前

https://www.onlinegdb.com/ に、突っ込んでみたところ、

コード:

#include <stdio.h>
#include <math.h>

double dob (double,double,double,double,double,double,double,double,double,double);

//PDゲイン
double kp = 60; //位置ゲイン 
double kv = 12; //速度ゲイン 
     
//運動モデル
const double dt = 0.01;    //サンプリングタイム
const double kt = 33.0;     //推力定数
const double omega = 500;   //カットオフ周波数 100

double tor_fb[3] = {0,0,0}; //DOB feedback
double tor4_old[3] = {0,0,0};
double tor5_old[3] = {0,0,0};


//エンドパートパラメータ
double L_AC = 0.07;
double L_BE = 0.1;
double L_HF = 0.1;
double L_HG = 0.07;
double L_OH = 0.13;

double MC = 0.05;
double ME = 0.07;
double MF = 0.05;
double MG = 0.07;
double MH = 0.07;





int main(void) {
    // file = fopen("/Users/tamuragou/Desktop/sample.txt","w");
    
    double t;                             //時間
    double g = 9.8;                       //重力加速度
    double pie = 3.1416;                  //π
    
    double xw_cmd, yw_cmd, zw_cmd;
    double dxw_cmd, dyw_cmd, dzw_cmd;;
    
    double ddxw_ref, ddyw_ref, ddzw_ref;
    
    double xw_res = 0, yw_res = -0.00707, zw_res = 0.137;
    double dxw_res = 0, dyw_res = 0, dzw_res =0;
    
    double j[3][3];
    double j_inv[3][3];
    
    double d = 0.01, the = pie/4, phi = 0;                      //モータ角度(初期値)
    double delta;
    
    double ddd_ref, ddthe_ref, ddphi_ref;
    
    double A,B;
    
    double m[3];
    double h_tor[3];
    double g_tor[3];
    
    double    d_res = 0,     the_res = 0,    phi_res = 0;
    double   dd_res = 0,    dthe_res = 0,   dphi_res = 0;
    double  ddd_res = 0,   ddthe_res = 0,  ddphi_res = 0;
    
    double    d_res_old = 0,     the_res_old = 0,    phi_res_old = 0;
    double   dd_res_old = 0,    dthe_res_old = 0,   dphi_res_old = 0;
    double  ddd_res_old = 0,   ddthe_res_old = 0,  ddphi_res_old = 0;
    
    double tor_1[3];
    double tor_2[3];
    double tor_3[3];
    
    double ia[3];
    
    double t_dis[3];
    double t_dis_hat[3];
    
    double ia_cmp[3] = {0,0,0};
    
    
    
    
    
    for (t = 0;t <= 2;t += dt){
        
        
        
        //Target position

            xw_cmd = 0.0001*t;
            yw_cmd = -0.01*sin(pie/4 + 0.01*t);
            zw_cmd = 0.13 + 0.01*cos(pie/4 + 0.01*t);

            dxw_cmd =  0.0001;
            dyw_cmd = -0.001*cos(0.785398 + 0.01*t);
            dzw_cmd = -0.001*sin(0.785398 + 0.01*t);

        
        
        
        //controller(PD制御)
        ddxw_ref = kp*(xw_cmd - xw_res)+kv*(dxw_cmd - dxw_res);
        ddyw_ref = kp*(yw_cmd - yw_res)+kv*(dyw_cmd - dyw_res);
        ddzw_ref = kp*(zw_cmd - zw_res)+kv*(dzw_cmd - dzw_res);
        
        
        //Forward Kinematics
        //ヤコビアンJ
        j[0][0] = sin(the)*sin(phi);
        j[0][1] = d*cos(the)*sin(phi);
        j[0][2] = d*sin(the)*cos(phi);
        
        j[1][0] = -sin(the)*cos(phi);
        j[1][1] = -d*cos(the)*cos(phi);
        j[1][2] =  d*sin(the)*sin(phi);
        
        j[2][0] = cos(the);
        j[2][1] = -d*sin(the);
        j[2][2] = 0;
        
        
        //Inverse Kinematics
        delta = j[0][1]*j[1][2]*j[2][0] + j[0][2]*j[1][0]*j[2][1] - j[0][2]*j[1][1]*j[2][0] - j[0][0]*j[1][2]*j[2][1];
        if( isnan( delta ) ) {
            printf( "delta is NaN %d.\n", __LINE__ );
            return -1;
        }
        
        if( isnan( -j[1][2] ) ) {
            printf( "-j[1][2] is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( j[2][1] ) ) {
            printf( "j[2][1] is NaN %d.\n", __LINE__ );
            return -1;
        }
        j_inv[0][0] = ( -j[1][2]*j[2][1] ) / delta;  // ここで、NaN 発生
        j_inv[0][1] = ( j[0][2]*j[2][1] )  / delta;
        j_inv[0][2] = ( j[0][1]*j[1][2] - j[0][2]*j[1][1] ) /delta;
        
        j_inv[1][0] =  ( j[1][2]*j[2][0] )  / delta;
        j_inv[1][1] =  ( -j[0][2]*j[2][0] ) / delta;
        j_inv[1][2] = -( j[0][0]*j[1][2] - j[0][2]*j[1][0] ) /delta;
        
        j_inv[2][0] =  ( j[1][0]*j[2][1] - j[1][1]*j[2][0] ) /delta;
        j_inv[2][1] = -( j[0][0]*j[2][1] - j[0][1]*j[2][0] ) /delta;
        j_inv[2][2] =  ( j[0][0]*j[1][1] - j[0][1]*j[1][0] ) /delta;
        
        
        //θ” = J^-1 * x"
        if( isnan( j_inv[0][0] ) ) {
            printf( "j_inv[0][0] is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( j_inv[0][1] ) ) {
            printf( "j_inv[0][1] is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( j_inv[0][2] ) ) {
            printf( "j_inv[0][2] is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( ddxw_ref ) ) {
            printf( "ddxw_ref is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( ddyw_ref ) ) {
            printf( "ddyw_ref is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( ddzw_ref ) ) {
            printf( "ddzw_ref is NaN %d.\n", __LINE__ );
            return -1;
        }
        ddd_ref   = j_inv[0][0]*ddxw_ref + j_inv[0][1]*ddyw_ref + j_inv[0][2]*ddzw_ref;
        ddthe_ref = j_inv[1][0]*ddxw_ref + j_inv[1][1]*ddyw_ref + j_inv[1][2]*ddzw_ref;
        ddphi_ref = j_inv[2][0]*ddxw_ref + j_inv[2][1]*ddyw_ref + j_inv[2][2]*ddzw_ref;
        
        
        A = MC*L_AC*L_AC + ME*L_BE*L_BE + MF*L_HF*L_HF + MG*L_HG*L_HG;
        B = MC*L_AC + ME*L_BE + MF*L_HF + MG*L_HG;
        
        
        //Inverse dynamics T = Mθ” + H + G
        //M
        m[0] = MH;
        m[1] = A + MH*d*d;
        m[2] = ( A + MH*d*d )*sin(the)*sin(the);
        
        //H
        h_tor[0] = -MH*d*( ( dthe_res*dthe_res ) + ( sin(the)*sin(the) )*( dphi_res*dphi_res ) );
        h_tor[1] = -( A + MH*d*d )*sin(the)*cos(the)*( dphi_res*dphi_res );
        h_tor[2] = 0;
        
        //G
        g_tor[0] = -MH*g*sin(the)*cos(phi);
        g_tor[1] = ( B - MH*d )*g*cos(the)*cos(phi);
        g_tor[2] = -(B - MH*d )*g*sin(the)*sin(phi);
        
        //T1 = M*θ”ref
        if( isnan( m[0] ) ) {
            printf( "m[0] is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( ddd_ref ) ) {
            printf( "ddd_ref is NaN %d.\n", __LINE__ );
            return -1;
        }
        tor_1[0] = m[0]*ddd_ref;
        tor_1[1] = m[1]*ddthe_ref;
        tor_1[2] = m[2]*ddphi_ref;
        
        //Ia = T1/Kt
        if( isnan( tor_1[0] ) ) {
            printf( "tor_1[0] is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( kt ) ) {
            printf( "kt is NaN %d.\n", __LINE__ );
            return -1;
        }
        ia[0] = tor_1[0]/kt;
        ia[1] = tor_1[1]/kt;
        ia[2] = tor_1[2]/kt;
        
        
        //外乱補償電流を加える
        if( isnan( ia[0] ) ) {
            printf( "ia[0] is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( ia_cmp[0] ) ) {
            printf( "ia_cmp[0] is NaN %d.\n", __LINE__ );
            return -1;
        }
        ia[0] += ia_cmp[0];
        ia[1] += ia_cmp[1];
        ia[2] += ia_cmp[2];
        
        
        //T2 = Kt*Ia
        if( isnan( kt ) ) {
            printf( "kt is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( ia[0] ) ) {
            printf( "ia[0] is NaN %d.\n", __LINE__ );
            return -1;
        }
        tor_2[0] = kt*ia[0];
        tor_2[1] = kt*ia[1];
        tor_2[2] = kt*ia[2];
        
        //外乱 Tdis = H + G
        t_dis[0] = h_tor[0] + g_tor[0];
        t_dis[1] = h_tor[1] + g_tor[1];
        t_dis[2] = h_tor[2] + g_tor[2];
        
        
        //T3 = T2 + Tdis
        if( isnan( tor_2[0] ) ) {
            printf( "tor_2[0] is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( t_dis[0] ) ) {
            printf( "t_dis[0] is NaN %d.\n", __LINE__ );
            return -1;
        }
        tor_3[0] = tor_2[0] + t_dis[0];
        tor_3[1] = tor_2[1] + t_dis[1];
        tor_3[2] = tor_2[2] + t_dis[2];
        
        
        //θ"res = T3/M
        if( isnan( tor_3[0] ) ) {
            printf( "tor_3 is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( m[0] ) ) {
            printf( "m[0] is NaN %d.\n", __LINE__ );
            return -1;
        }
        ddd_res   = tor_3[0]/m[0];
        if( isnan( ddd_res ) ) {
            printf( "ddd_res is NaN %d.\n", __LINE__ );
            return -1;
        }
        ddthe_res = tor_3[1]/m[1];
        ddphi_res = tor_3[2]/m[2];
        
        //θres = 1/s * θ'res
        if( isnan( dd_res ) ) {
            printf( "dd_res is NaN %d.\n", __LINE__ );
            return -1;
        }
        d_res = d_res_old   + (dd_res   + dd_res_old)*dt/2.0;
        if( isnan( d_res ) ) {
            printf( "d_res is NaN %d.\n", __LINE__ );
            return -1;
        }

        the_res = the_res_old + (dthe_res + dthe_res_old)*dt/2.0;
        phi_res = phi_res_old + (dphi_res + dphi_res_old)*dt/2.0;
        
        
        //θ'res = 1/s * θ"res
        if( isnan( dd_res_old ) ) {
            printf( "dd_res is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( ddd_res ) ) {
            printf( "ddd_res is NaN %d.\n", __LINE__ );
            return -1;
        }
        if( isnan( ddd_res_old ) ) {
            printf( "ddd_res_old is NaN %d.\n", __LINE__ );
            return -1;
        }
        dd_res   = dd_res_old   + (ddd_res   + ddd_res_old)*dt/2.0;
        if( isnan( dd_res ) ) {
            printf( "dd_res is NaN %d.\n", __LINE__ );
            return -1;
        }
        dthe_res = dthe_res_old + (ddthe_res + ddthe_res_old)*dt/2.0;
        dphi_res = dphi_res_old + (ddphi_res + ddphi_res_old)*dt/2.0;
       
        
        //前回値としてθ"resを記憶
        ddd_res_old   = ddd_res;
        ddthe_res_old = ddthe_res;
        ddphi_res_old = ddphi_res;
        
        //前回値としてθ'resを記憶
        dd_res_old   = dd_res;
        dthe_res_old = dthe_res;
        dphi_res_old = dphi_res;
        
        //前回値としてθresを記憶
        d_res_old   = d_res;
        the_res_old = the_res;
        phi_res_old = phi_res;
        
        
        //Tdis^ (DOB)
        t_dis_hat[0] = dob(ia[0], ia[1], ia[2], dd_res, dthe_res, dphi_res, m[0], m[1], m[2], 0);
        t_dis_hat[1] = dob(ia[0], ia[1], ia[2], dd_res, dthe_res, dphi_res, m[0], m[1], m[2], 1);
        t_dis_hat[2] = dob(ia[0], ia[1], ia[2], dd_res, dthe_res, dphi_res, m[0], m[1], m[2], 2);
        
        
        //外乱補償電流Iacmp
        ia_cmp[0] = (1/kt) * t_dis_hat[0];
        ia_cmp[1] = (1/kt) * t_dis_hat[1];
        ia_cmp[2] = (1/kt) * t_dis_hat[2];
        
        
        //x'res = J * θ'res
        dxw_res = j[0][0]*dd_res + j[0][1]*dthe_res + j[0][2]*dphi_res;
        dyw_res = j[1][0]*dd_res + j[1][1]*dthe_res + j[1][2]*dphi_res;
        dzw_res = j[2][0]*dd_res + j[2][1]*dthe_res + j[2][2]*dphi_res;
        
        
        //xresをθresから求める
        xw_res =  d_res*sin(the_res)*sin(phi_res);
        yw_res = -d_res*sin(the_res)*cos(phi_res);
        
        if( isnan( d_res ) ) {
            printf( "d_res is NaN %d.\n", __LINE__ );
            return -1;
        }
        zw_res = L_OH + d_res*cos(the_res);
        
        
        d = d_res;
        the = the_res;
        phi = phi_res;

        printf("%lf\t%lf \n",t,zw_res);
        //fprintf(file,"%lf\t%lf  \n",t,zw_res);
    }
    
}


//外乱オブザーバー
double dob (double ia_dob0, double ia_dob1, double ia_dob2, double d_dob, double the_dob, double phi_dob, double m0, double m1, double m2, double count){
    
    double tor_d1[3];
    double tor_d2[3];
    double tor_d3[3];
    double tor_d4[3];
    double tor_d5[3];
    double tdis_dob[3];
    
    
    
    //T1 = kt*Ia
    tor_d1[0] = kt*ia_dob0;
    tor_d1[1] = kt*ia_dob1;
    tor_d1[2] = kt*ia_dob2;
    
    
    //T2 = ωc*(MN * θ'res)
    tor_d2[0] = omega*( m0*d_dob );
    tor_d2[1] = omega*( m1*the_dob );
    tor_d2[2] = omega*( m2*phi_dob );
    
    
    //T3 = T1 + T2
    tor_d3[0] = tor_d1[0] + tor_d2[0];
    tor_d3[1] = tor_d1[1] + tor_d2[1];
    tor_d3[2] = tor_d1[2] + tor_d2[2];
    
    
    //T4 = ωc*T3 - T_feedback
    tor_d4[0] = ( omega*tor_d3[0] ) - tor_fb[0];
    tor_d4[1] = ( omega*tor_d3[1] ) - tor_fb[1];
    tor_d4[2] = ( omega*tor_d3[2] ) - tor_fb[2];
    
    
    //LPF///////////////////////////////////////
    //T5 = 1/s * T4
    tor_d5[0] = tor5_old[0] + ( tor_d4[0] + tor4_old[0] )*dt/2.0;
    tor_d5[1] = tor5_old[1] + ( tor_d4[1] + tor4_old[1] )*dt/2.0;
    tor_d5[2] = tor5_old[2] + ( tor_d4[2] + tor4_old[2] )*dt/2.0;
    
    
    //T_feedback = ωc * T5 (フィードバック)
    tor_fb[0] = omega*tor_d5[0];
    tor_fb[1] = omega*tor_d5[1];
    tor_fb[2] = omega*tor_d5[2];
    
    
    //前回値としてT4を記憶
    tor4_old[0] = tor_d4[0];
    tor4_old[1] = tor_d4[1];
    tor4_old[2] = tor_d4[2];
    
    
    //前回値としてT5を記憶
    tor5_old[0] = tor_d5[0];
    tor5_old[1] = tor_d5[1];
    tor5_old[2] = tor_d5[2];
    
    //.//////////////////////////////////////////
    
    
    //Tdis = T5 - T2
    tdis_dob[0] =  tor_d5[0] - tor_d2[0];
    tdis_dob[1] =  tor_d5[1] - tor_d2[1];
    tdis_dob[2] =  tor_d5[2] - tor_d2[2];
    
    
    //出力
       if(count == 0){
           return tdis_dob[0];
           
       }else if (count == 1){
           return tdis_dob[1];
           
       }else if (count == 2){                     //default
           return tdis_dob[2];
           
       }else{                 //default
           return tdis_dob[2];
       }
    
}
結果

コード:

0.000000         0.130000
J_inv[0][0] is NaN 157.
142 行目の
j_inv[0][0] = ( -j[1][2]*j[2][1] ) / delta; // ここで、NaN 発生
VTuber:
東上☆海美☆(とうじょう・うみみ)
http://atassyu.php.xdomain.jp/vtuber/index.html
レスがついていないものを優先して、レスするみみ。時々、見当外れなレスしみみ。

中の人:
手提鞄あたッしュ、[MrAtassyu] 手提鞄屋魚有店
http://ameblo.jp/mratassyu/
Pixiv: 666303
Windows, Mac, Linux, Haiku, Raspbery Pi, Jetson Nano, 電子ブロック 持ち。

アバター
lego hasiri
記事: 7
登録日時: 6年前
連絡を取る:

Re: NaNが出る原因が何か分からない

#3

投稿記事 by lego hasiri » 3年前

簡単に見ただけですけど、forの2週目でdとかtheとかの値がすべて0になっていると思います。1週目は動いていそうなので、for後半のパラメータの計算過程が間違っていそうかな?

返信

“C言語何でも質問掲示板” へ戻る