Eigen初心者ですが,状態空間方程式を解く際にC++でEigenを用いてプログラムを作成しているのですが,なかなかうまくいきません.解法も含めて教えていただきたいです.
問題は以下の通りですが,
- この系の数理模型を作成せよ。
2. この系の状態空間モデルを作成せよ。ただし、状態x(t) は
以下のように定義せよ。
x(t) = col(y(t), y˙(t))
3. この系のデジタルシミュレータを作成せよ。
4. この系をevent(t0, x0) からevent(t1, x1) に移すための制御入力を見出せ。
ただし、時刻t0 = 0、状態x0 = col(1, 0)、時刻t1 = 1、
状態x1 = col(0, 0) とせよ。
自分での解釈は以下の通りなのですが,
これが間違ってたら指摘お願いします.
ソースコードは以下の通りです.
なかなかうまくいかなくて困ってます.教えていただけるとありがたいです.
C++
1#include "Eigen/Dense" 2#include "Eigen/Core" 3#include <iostream> 4#include <fstream> 5#include <cmath> 6using namespace Eigen; 7using namespace std; 8 9void RungeKutta(MatrixXd &X, double tt, double dt, MatrixXd A, MatrixXd B, MatrixXd U) { 10 MatrixXd k1 = A*X + B*U; 11 MatrixXd k2 = A*(X + 0.5*k1*dt) + B*U; 12 MatrixXd k3 = A*(X + 0.5*k2*dt) + B*U; 13 MatrixXd k4 = A*(X + k3*dt) + B*U; 14 MatrixXd k = (k1 + 2.0*k2 + 2.0*k3 + k4)*dt / 6.0; 15 X = X + k; 16} 17 18int main() { 19 double m = 1.0; 20 21 MatrixXd A(2, 2); 22 A(0, 0) = 0; 23 A(0, 1) = 0; 24 A(1, 0) = 0; 25 A(1, 1) = 1; 26 MatrixXd B(2, 2); 27 B(0, 0) = 0; 28 B(0, 1) = 0; 29 B(1, 0) = 0; 30 B(1, 1) = 1 / m; 31 32 double dt = 0.01; 33 double tt = 0.0; 34 MatrixXd X(2, 1); 35 X(0, 0) = 1; 36 X(1, 0) = 0; 37 MatrixXd U(2, 1); 38 U(0, 0) = 0; 39 U(1, 0) = -6.25*sin(2.0*tt*M_PI); 40 41 double freq = 1.0; 42 ofstream ofs("outdata.csv"); 43 ofs << "time," << "y" << endl; 44 45 for (int i = 0; i < 1000; i++) { 46 RungeKutta(X, tt, dt, A, B, U); 47 ofs << tt << "," << X(0, 0) << endl; 48 tt += dt; 49 U(0, 1) = -6.25*sin(2.0*tt*M_PI); 50 } 51 52 return 0; 53}
エラー内容は,
Assertion failed: row >= 0 && row < rows() && col >= 0 && col < cols(), file Eigen/src/Core/DenseCoeffsBase.h, line 365
です.
よろしくお願いいたします.
追記!
c++
1#include "Eigen/Dense" 2#include "Eigen/Core" 3#include <iostream> 4#include <fstream> 5#include <cmath> 6using namespace Eigen; 7using namespace std; 8 9 10void RungeKutta(MatrixXd &X, double tt, double dt, MatrixXd A, MatrixXd B, MatrixXd U) { 11 MatrixXd k1 = A*X*dt + -6.25*sin(2.0*tt*M_PI)*B*U*dt; 12 MatrixXd k2 = A*(X + 0.5*k1)*dt + -6.25*sin(2.0*tt*M_PI)*B*U*dt; 13 MatrixXd k3 = A*(X + 0.5*k2)*dt + -6.25*sin(2.0*tt*M_PI)*B*U*dt; 14 MatrixXd k4 = A*(X + k3)*dt + -6.25*sin(2.0*tt*M_PI)*B*U*dt; 15 MatrixXd k = (k1 + 2.0*k2 + 2.0*k3 + k4)*dt / 6.0; 16 X = X + k; 17} 18 19int main() { 20 double m = 1.0; 21 22 MatrixXd A(2, 2); 23 A(0, 0) = 0; 24 A(0, 1) = 1; 25 A(1, 0) = 0; 26 A(1, 1) = 0; 27 MatrixXd B(2, 2); 28 B(0, 0) = 0; 29 B(0, 1) = 0; 30 B(1, 0) = 0; 31 B(1, 1) = 1 / m; 32 33 double dt = 0.001; 34 double tt = 0.0; 35 MatrixXd X(1, 2); 36 X(0, 0) = 1; 37 X(0, 1) = 0; 38 MatrixXd U(1, 2); 39 U(0, 0) = 0; 40 U(0, 1) = 1; 41 42 double freq = 1.0; 43 ofstream ofs("outdata.csv"); 44 ofs << "time," << "x," << "v" << endl; 45 46 for (int i = 0; i <= 1000; i++) { 47 RungeKutta(X, tt, dt, A, B, U); 48 ofs << tt << "," << X(0, 0) << "," << X(0, 1) << endl; 49 tt += dt; 50 } 51 52 return 0; 53}
回答1件
あなたの回答
tips
プレビュー
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。
2019/01/07 12:11
2019/01/07 12:16
2019/01/07 12:17
2019/01/07 12:18
2019/01/07 12:20
2019/01/07 12:37
2019/01/07 12:38
2019/01/08 08:04
2019/01/09 23:54