制振時と非制振時における波形計算で，同調調和減衰器（TMD）の制振効果を比較するMaTXプログラム．
なお，本プログラムは，MaTX公式サイトのMaTXユーザーズマニュアルの例題プログラムを参考にしています．

グラフは，gnuplotにより以下のように出力されます．

// feedback.mm

Matrix Ap, bp, cp;     // 制御対象のパラメータ
Matrix Ac, bc, cc;     // 制御器のパラメータ
Matrix As, bs;         // 制御対象（非制御時）のパラメータ

Func void init_params()
{
  Real f_p, d_p, m_p, f_c, d_c, m_c;

  f_p=2*PI/2.70;
  d_p=0.0005;
  m_p=1200;
  f_c=2*PI/2.72;
  d_c=0.0660;
  m_c=15.4;

  // 制御対象のパラメータ
  Ap = [[0, 1]
        [-f_p^2-(m_c/m_p)*f_c^2, -2*d_p*f_p-2*(m_c/m_p)*d_c*f_c]];
  bp = [0 1/m_p]';
  cp = [m_c*f_c^2, 2*m_c*d_c*f_c];
  // 制御器のパラメータ
  Ac = [[0,  1]
        [-f_c^2,  -2*d_c*f_c]];
  bc = [0 1/m_c]';
  cc = [m_c*f_c^2, 2*m_c*d_c*f_c];
  // 制御対象のパラメータ（非制御時）
  As = [[0, 1]
        [-f_p^2, -2*d_p*f_p]];
  bs = [0 1/m_p]';
}

Func void main()
{
  Matrix x0;
  Matrix xs0;
  Array TC, XC, UC;
   
  void diff_eqs(), link_eqs(), init_params();

  init_params();

  x0 = [0 0 0 0 0 0]';
  {TC, XC, UC} =
     Ode(0.0, 100.0, x0, diff_eqs, link_eqs, 0.1);
  mgplot(1, TC, [[XC(1:1,:)][XC(3:3,:)][XC(5:5,:)]], {"xp_disp", "xc_disp", "xs_disp"});
}

Func void diff_eqs(DX, t, X, UY)
  Real t;
  Matrix X, DX, UY;
{
  Matrix dxp, dxc, dxs, xp, xc, xs, up, uc, us;

  xp = X(1:2,1);          // 制御対象の状態
  xc = X(3:4,1);          // 制御器の状態
  xs = X(5:6,1);          // 制御対象の状態（非制御時）
  up = UY(1:1,1);         // 制御対象への入力
  uc = UY(2:2,1);         // 制御器への入力
  us = UY(3:3,1);         // 制御対象への入力（非制御時）

  dxp = Ap*xp + bp*up;    // 制御対象の状態方程式
  dxc = Ac*xc + bc*uc;    // 制御器の状態方程式
  dxs = As*xs + bs*us;    // 制御対象の状態方程式（非制御時）

  DX = [[dxp][dxc][dxs]];
}

Func void link_eqs(UY, t, X)
  Real t;
  Matrix UY, X;
{
  Matrix xp, up, yp, xc, yc, uc, xs, us;
  Real noise;

  xp = X(1:2,1);          // 制御対象の状態
  xc = X(3:4,1);          // 制御器の状態
  xs = X(5:6,1);          // 制御対象の状態（非制御時）

  yp = cp*xp;             // 制御対象の出力
  uc = yp;                // yp を uc に接続

  yc = cc*xc;             // 制御器の出力
  us = [(rand() - 0.5)*2];
  up = us + yc;           // yc を up に接続

  UY = [up uc us]';
}