MATLAB,Python,Scilab,Julia比較 その62【状態空間モデル⑳】

バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia-backnumber/

はじめに
登場人物
Cコード化可能か？
Cコード
シミュレーション結果
まとめ

はじめに

前回までで、各ツール各言語による状態空間モデルのベクトル行列演算による実現をしたところ。

よくよく考えると、Cコード化という作業があるようなないような。

登場人物

博識フクロウのフクさん

イラストACにて公開の「kino_k」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=iKciwKA9&area=1

エンジニア歴8年の太郎くん

イラストACにて公開の「しのみ」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=uCKphAW2&area=1

Cコード化可能か？

太郎くん

ふと思ったんだけど？

フクさん

何？

太郎くん

状態空間モデルを使用しない場合はCコード化までやったじゃん？
状態空間モデルを使用する場合Cコード化はどうなるの？

フクさん

まぁベクトル行列演算用のライブラリを作ってやることもあるが、
C言語の場合、動的に配列を定義することが難しかったり、
オーバーヘッドがあったりする。
よって、想定される配列の大きさを決めてからCコード化することは多いな。

太郎くん

うーん、イメージがわかん。

フクさん

じゃ、とりあえずMATLAB Coderで今回の状態空間モデルがどのようにCコードとして表現されるか見ているか。

太郎くん

自分で書かないのか？！

フクさん

いや、結構メンドクサイから・・・。

フクさん

まぁ、まずは微分解決済みの状態空間モデルの再掲しておこう。

状態方程式
\(\boldsymbol{x}(t+\Delta t)=\boldsymbol{x}(t)+\{A\boldsymbol{x}(t)+B\boldsymbol{u}(t)\}\Delta t \)

出力方程式
\(\boldsymbol{y}(t+\Delta T)=C\boldsymbol{x}(t+\Delta t)+D\boldsymbol{u}(t)\)

Cコード

フクさん

以下がMATLAB Coderが出力したものにたいして、若干の修正を入れたものだ。

/*
 * main.c
 *
 * Code generation for function 'main'
 *
 */

/* Include files */
#include <stdio.h>
#include "rt_nonfinite.h"
#include "main.h"

/* Type Definitions */

/* Named Constants */

/* Variable Declarations */

/* Variable Definitions */

/* Function Declarations */
static void b_statespacemodel(const real_T A[4], const real_T B[2], const real_T
  C[4], const real_T D[2], real_T u, real_T dt, real_T x[2], real_T y[2]);
static void eml_xgemm(const real_T A[4], const real_T B[2], real_T C[2]);

/* Function Definitions */
static void eml_xgemm(const real_T A[4], const real_T B[2], real_T C[2])
{
  int32_T i1;
  int32_T i2;
  for (i1 = 0; i1 < 2; i1++) {
    C[i1] = 0.0;
    for (i2 = 0; i2 < 2; i2++) {
      C[i1] += A[i1 + (i2 << 1)] * B[i2];
    }
  }
}

real_T y[20000];

void main()
{
  statespacemodel_motion_raw(y);
}

void statespacemodel(const real_T A[4], const real_T B[2], const real_T C[4],
                     const real_T D[2], real_T u, real_T dt, real_T x[2], real_T
                     y[2])
{
  real_T dv5[2];
  int32_T i;

  /*  様態方程式 */
  for (i = 0; i < 2; i++) {
    dv5[i] = 0.0;
  }

  eml_xgemm(A, x, dv5);
  for (i = 0; i < 2; i++) {
    x[i] += (dv5[i] + B[i] * u) * dt;
  }

  /*  出力方程式 */
  for (i = 0; i < 2; i++) {
    y[i] = 0.0;
  }

  eml_xgemm(C, x, y);
  for (i = 0; i < 2; i++) {
    y[i] += D[i] * u;
  }
}

int8_T u[10000];

void statespacemodel_motion_raw()
{
  real_T x[2];
  int32_T i;
  real_T dv0[2];
  real_T dv1[2];
  int32_T i0;
  static const real_T dv2[4] = { 1.0, 0.0, 0.0, 1.0 };

  static const real_T dv3[4] = { 0.0, 1.0, 0.0, 0.0 };

  /*  時間(横)軸 */
  memset(&u[0], 0, 10000U * sizeof(int8_T));

  /*  入力信号生成 */
  memset(&u[4999], 1, 5001U * sizeof(int8_T));

  /*  5秒後に0から1へ */
  memset(&y[0], 0, 20000U * sizeof(real_T));
  for (i = 0; i < 2; i++) {
    x[i] = 0.0;
  }
  
  printf("F,v,s\n");
  for (i = 0; i < 10000; i++) {
    for (i0 = 0; i0 < 2; i0++) {
      dv0[i0] = 1.0 - (real_T)i0;
      dv1[i0] = 0.0;
    }

    statespacemodel(dv3, dv0, dv2, dv1, (real_T)u[i], 0.001, x, *(real_T (*)[2])
                      &y[i << 1]);
    printf("%d,%lf,%lf\n",u[i], x[0], x[1]);
  }
}

/* End of code generation (main.c) */

太郎くん

修正ってなにしたの？

フクさん

ローカル変数に巨大な配列が定義されてしまうので、
それをグローバル変数に移動させた。
以下が該当する。