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

グラフは，gnuplotにより以下のように出力されます．緑線が制振なしの場合の変位波形，青線が制振ありの場合の変位波形です．

// tmd.mm

Matrix Ap, bp, cp;     // 制御対象のパラメータ
Matrix As, bs, cs;     // 制御対象（非制御時）のパラメータ
Real f_s, d_s, m_s, f_d, d_d, m_d;

Func void init_params()
{
  Matrix MM, CC, KK;

  f_s=2*PI/2.70;
  d_s=0.0005;
  m_s=1200;
  f_d=2*PI/2.72;
  d_d=0.0660;
  m_d=15.4;
  
  MM = [[m_s, 0]
        [0, m_d]];
  CC = [[2*m_s*d_s*f_s+2*m_d*d_d*f_d, -2*m_d*d_d*f_d]
        [-2*m_d*d_d*f_d, 2*m_d*d_d*f_d]];
  KK = [[m_s*f_s^2+m_d*f_d^2, -m_d*f_d^2]
        [-m_d*f_d^2, m_d*f_d^2]];

  // 制御対象のパラメータ
  Ap = [[Z(2), diag(1.0, 1.0)]
        [-MM~*KK, -MM~*CC]];
  bp = [[Z(2)]
        [MM~]];
  cp = [[0, 0, 0, 0]
        [0, 0, 0, 0]];
  // 制御対象のパラメータ（非制御時）
  As = [[0, 1]
        [-f_s^2, -2*d_s*f_s]];
  bs = [0, 1/m_s]';
  cs = [0, 0];
}

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

  init_params();

  print Ap;
  print eigval(Ap);
  print As;
  print eigval(As);

  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(5:5,:)]], {"with TMD", "without TMD"});
  pause;
}

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

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

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

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

Func void link_eqs(UY, t, X)
  Real t;
  Matrix UY, X;
{
  Matrix xp, up, yp, xs, us, ys, vp, vs;
  Real noise;

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

  noise = (rand() - 0.5)*2;

  vp = [noise 0]';        // 外乱
  vs = [noise];
  yp = cp*xp;             // 出力
  ys = cs*xs;
  up = yp+vp;             // 入力
  us = ys+vs;

  UY = [[up][us]];
}