Dlsim.java

/*
 * $Id: Dlsim.java,v 1.24 2008/03/23 14:29:07 koga Exp $
 *
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.control;

import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.tool.graph.gnuplot.Canvas;
import org.mklab.tool.graph.gnuplot.Gnuplot;


/**
 * 離散時間線形システム任意入力に対する応答を求めるクラスです。
 * 
 * <p>Response of disc-time linear systems to given input
 * 
 * @author koga
 * @version $Revision: 1.24 $
 * @see org.mklab.tool.control.Lsim
 */
public class Dlsim {

  /**
   * 離散時間システム
   * 
   * <pre><code>
   * 
   * x[n+1] = Ax[n] + Bu[n] y[n] = Cx[n] + Du[n]
   * 
   * </code></pre>
   * 
   * に入力列 <code>U</code> を与え、その時間応答を計算します。 行列 <code>U</code> の各列が各時刻の入力です。
   * 
   * @param A A行列
   * @param B B行列
   * @param C C行列
   * @param D D行列
   * @param U 入力の列
   * @return {Y, x} (出力と状態応答) response
   */
  public static List<DoubleMatrix> dlsim(DoubleMatrix A, DoubleMatrix B, DoubleMatrix C, DoubleMatrix D, DoubleMatrix U) {
    DoubleMatrix x0 = A.createZero(A.getRowSize(), 1);
    return dlsim(A, B, C, D, U, x0);
  }

  /**
   * <code>x0</code>を初期状態として指定します。
   * 
   * @param A A行列
   * @param B B行列
   * @param C C行列
   * @param D D行列
   * @param U 入力の列
   * @param x0 初期状態
   * @return {Y, x} (出力と状態応答) response
   */
  public static List<DoubleMatrix> dlsim(DoubleMatrix A, DoubleMatrix B, DoubleMatrix C, DoubleMatrix D, DoubleMatrix U, DoubleMatrix x0) {
    String message;
    if ((message = Abcdchk.abcdchk(A, B, C, D)).length() > 0) {
      throw new RuntimeException(message);
    }

    DoubleMatrix X = Ltitr.ltitr(A, B, U, x0);
    DoubleMatrix Y = C.multiply(X).add(D.multiply(U));
    return new ArrayList<>(Arrays.asList(new DoubleMatrix[] {Y, X}));
    //return new MatxList(new Object[] {Y, X});
  }

  /**
   * @param A A行列
   * @param B B行列
   * @param C C行列
   * @param D D行列
   * @param U 入力の列
   * @return Gnuplot
   * @throws IOException gnuplotプロセスを起動できない場合
   */
  public static Gnuplot plot(DoubleMatrix A, DoubleMatrix B, DoubleMatrix C, DoubleMatrix D, DoubleMatrix U) throws IOException {
    Gnuplot gp = new Gnuplot();
    return plot(gp, A, B, C, D, U);
  }

  /**
   * @param gnuplot gnuplot
   * @param A A行列
   * @param B B行列
   * @param C C行列
   * @param D D行列
   * @param U 入力の列
   * @return Gnuplot
   */
  public static Gnuplot plot(Gnuplot gnuplot, DoubleMatrix A, DoubleMatrix B, DoubleMatrix C, DoubleMatrix D, DoubleMatrix U) {
    DoubleMatrix x0 = new DoubleMatrix(A.getRowSize(), 1);
    return plot(gnuplot, A, B, C, D, U, x0);
  }

  /**
   * @param A A行列
   * @param B B行列
   * @param C C行列
   * @param D D行列
   * @param U 入力の列
   * @param x0 初期状態
   * @return Gnuplot
   * @throws IOException gnuplotプロセスを起動できない場合
   */
  public static Gnuplot plot(DoubleMatrix A, DoubleMatrix B, DoubleMatrix C, DoubleMatrix D, DoubleMatrix U, DoubleMatrix x0) throws IOException {
    Gnuplot gp = new Gnuplot();
    return plot(gp, A, B, C, D, U, x0);
  }

  /**
   * @param gnuplot gnuplot
   * @param A A行列
   * @param B B行列
   * @param C C行列
   * @param D D行列
   * @param U 入力の列
   * @param x0 初期状態
   * @return Gnuplot
   */
  @SuppressWarnings("nls")
  public static Gnuplot plot(Gnuplot gnuplot, DoubleMatrix A, DoubleMatrix B, DoubleMatrix C, DoubleMatrix D, DoubleMatrix U, DoubleMatrix x0) {
    String msg = Abcdchk.abcdchk(A, B, C, D);
    if (msg.length() > 0) {
      throw new RuntimeException(msg);
    }

    List<DoubleMatrix> yx = dlsim(A, B, C, D, U, x0);
    DoubleMatrix Y = yx.get(0);
    DoubleMatrix X = yx.get(1);

    int outputSize = C.getRowSize();
    String[] yy = new String[outputSize];
    for (int i = 1; i <= outputSize; i++) {
      if (outputSize == 1) {
        yy[i - 1] = "y";
      } else {
        yy[i - 1] = "y" + i;
      }
    }

    gnuplot.createCanvas(2, 1);
    Canvas canvas1 = gnuplot.getCanvas(0, 0);
    canvas1.setGridVisible(true);
    canvas1.setXLabel("k");
    canvas1.setYLabel("y");

    int N = U.getColumnSize();

    canvas1.plot(DoubleMatrix.series(1, N, 1), Y, yy);

    int stateSize = A.getRowSize();
    String[] xx = new String[stateSize];
    for (int i = 1; i <= stateSize; i++) {
      xx[i - 1] = "x" + i;
    }

    Canvas canvas2 = gnuplot.getCanvas(1, 0);
    canvas2.setGridVisible(true);
    canvas2.setXLabel("k");
    canvas2.setYLabel("x");
    canvas2.plot(DoubleMatrix.series(1, N, 1), X, xx);
    return gnuplot;
  }

}