Impulse.java

/*
 * $Id: Impulse.java,v 1.31 2008/03/24 23:45:43 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.nfc.matrix.DoubleRationalPolynomialMatrix;
import org.mklab.nfc.scalar.DoubleNumber;
import org.mklab.nfc.scalar.DoubleRationalPolynomial;
import org.mklab.tool.graph.gnuplot.Canvas;
import org.mklab.tool.graph.gnuplot.Gnuplot;


/**
 * 連続系のインパルス応答を求めるクラスです。
 * 
 * <p>Impulse response of continuous-time linear system
 * 
 * @author koga
 * @version $Revision: 1.31 $
 * @see org.mklab.tool.control.Step
 * @see org.mklab.tool.control.Dimpulse
 */
public class Impulse {

  /**
   * 連続時間システム
   * 
   * <pre><code>
   * 
   * . x = Ax + Bu y = Cx + Du
   * 
   * </code></pre>
   * 
   * の <code>iu</code> 番目の入力に単位インパルスが加えられた場合の インパルス応答を計算します。
   * 
   * <p> <code>iu == 0</code> の場合、
   * 
   * <pre><code>
   * 
   * [[Y for 1st input] [[X for 1st input] [Y for 2nd input] [X for 2nd input] [...............] [...............] [Y for m'th input]] と [Y for m'th input]].
   * 
   * </code></pre>
   * 
   * を返します。
   * 
   * @param A A行列
   * @param B B行列
   * @param C C行列
   * @param D D行列
   * @param inputNumber 入力番号
   * @param T 評価する時刻の列(等間隔)
   * @return {YY, XX} (出力と状態の応答) response
   */
  public static List<DoubleMatrix> impulse(DoubleMatrix A, DoubleMatrix B, DoubleMatrix C, DoubleMatrix D, int inputNumber, DoubleMatrix T) {
    String message;
    if ((message = Abcdchk.abcdchk(A, B, C, D)).length() > 0) {
      throw new RuntimeException(message);
    }

    int N = T.length();
    DoubleMatrix U = A.createZero(1, N);
    DoubleNumber dt = T.getElement(2).subtract(T.getElement(1));

    DoubleMatrix XX, YY;

    if (inputNumber == 0) {
      int inputSize = B.getColumnSize();
      YY = A.createZero(C.getRowSize() * inputSize, N);
      XX = A.createZero(A.getRowSize() * inputSize, N);

      for (int i = 1; i <= inputSize; i++) {
        DoubleMatrix b = B.getColumnVector(i);
        DoubleMatrix x0 = b;
        List<DoubleMatrix> ab = C2d.c2d(A, b, dt);
        DoubleMatrix Ad = ab.get(0);
        DoubleMatrix Bd = ab.get(1);
        DoubleMatrix X = Ltitr.ltitr(Ad, Bd, U, x0);
        DoubleMatrix Y = C.multiply(X);

        YY.setSubMatrix(i, 1, Y, Y);
        XX.setSubMatrix(i, 1, X, X);
      }
    } else {
      DoubleMatrix b = B.getColumnVector(inputNumber);
      DoubleMatrix x0 = b;
      List<DoubleMatrix> ab = C2d.c2d(A, b, dt);
      DoubleMatrix Ad = ab.get(0);
      DoubleMatrix Bd = ab.get(1);
      XX = Ltitr.ltitr(Ad, Bd, U, x0);
      YY = C.multiply(XX);
    }

    return new ArrayList<>(Arrays.asList(new DoubleMatrix[] {YY, XX}));
    //return new MatxList(new Object[] {YY, XX});
  }

  /**
   * 単位インパルス応答を計算する
   * 
   * @param system 線形システム
   * @param T 評価する時刻の列
   * @return {YY, XX} (出力と状態の応答) response
   */
  public static List<DoubleMatrix> impulse(DoubleLinearSystem system, DoubleMatrix T) {
    if (system.isProper()) {
      DoubleMatrix A = system.getA();
      DoubleMatrix B = system.getB();
      DoubleMatrix C = system.getC();
      DoubleMatrix D = system.getD();
      return impulse(A, B, C, D, 1, T);
    }

    return impulse(system.getTransferFunctionMatrix(), 1, T);
  }

  /**
   * 単位インパルス応答を計算する
   * 
   * @param g 伝達関数
   * @param T 評価する時刻の列
   * @return {YY, XX} (出力と状態の応答) response
   */
  public static List<DoubleMatrix> impulse(DoubleRationalPolynomial g, DoubleMatrix T) {
    List<DoubleMatrix> abcd = Tfn2ss.tfn2ss(g);
    DoubleMatrix A = abcd.get(0);
    DoubleMatrix B = abcd.get(1);
    DoubleMatrix C = abcd.get(2);
    DoubleMatrix D = abcd.get(3);
    return impulse(A, B, C, D, 1, T);
  }

  /**
   * 単位インパルス応答を計算する
   * 
   * @param G 伝達関数行列
   * @param inputNumber 入力番号
   * @param T 評価する時刻の列
   * @return {YY, XX} (出力と状態の応答) response
   */
  public static List<DoubleMatrix> impulse(DoubleRationalPolynomialMatrix G, int inputNumber, DoubleMatrix T) {
    if (inputNumber == 0) {
      List<DoubleMatrix> abcd = Tfm2ss.tfm2ss(G);
      DoubleMatrix A = abcd.get(0);
      DoubleMatrix B = abcd.get(1);
      DoubleMatrix C = abcd.get(2);
      DoubleMatrix D = abcd.get(3);
      return impulse(A, B, C, D, inputNumber, T);
    }
    List<DoubleMatrix> abcd = Tfm2ss.tfm2ss(G, inputNumber);
    DoubleMatrix A = abcd.get(0);
    DoubleMatrix B = abcd.get(1);
    DoubleMatrix C = abcd.get(2);
    DoubleMatrix D = abcd.get(3);
    return impulse(A, B, C, D, 1, T);
  }

  /**
   * 単位インパルス応答を計算する
   * 
   * @param numerator 伝達関数の分子多項式の係数
   * @param denominator 伝達関数の分母多項式の係数
   * @param T 評価する時刻の列
   * @return {YY, XX} (出力と状態の応答) response
   */
  public static List<DoubleMatrix> impulse(DoubleMatrix numerator, DoubleMatrix denominator, DoubleMatrix T) {
    List<DoubleMatrix> abcd = Tf2ss.tf2ss(numerator, denominator);
    DoubleMatrix A = abcd.get(0);
    DoubleMatrix B = abcd.get(1);
    DoubleMatrix C = abcd.get(2);
    DoubleMatrix D = abcd.get(3);
    return impulse(A, B, C, D, 1, T);
  }

  /**
   * @param g 伝達関数
   * @param T 評価する時刻の列
   * @return Gnuplot
   * @throws IOException gnuplotプロセスを起動できない場合
   */
  public static Gnuplot plot(DoubleRationalPolynomial g, DoubleMatrix T) throws IOException {
    Gnuplot gp = new Gnuplot();
    return plot(gp, g, T);
  }

  /**
   * 単位インパルス応答をプロットする
   * 
   * @param gnuplot gnuplot
   * @param g 伝達関数
   * @param T 評価する時刻の列
   * @return Gnuplot
   */
  @SuppressWarnings("nls")
  public static Gnuplot plot(Gnuplot gnuplot, DoubleRationalPolynomial g, DoubleMatrix T) {
    List<DoubleMatrix> yx = impulse(g, T);
    DoubleMatrix Y = yx.get(0);
    Canvas canvas = gnuplot.createCanvas();
    canvas.setGridVisible(true);
    canvas.setXLabel("t [sec]");
    canvas.setYLabel("y");
    canvas.plot(T, Y, new String[] {"y"});
    return gnuplot;
  }

  /**
   * @param numerator 分子多項式の係数
   * @param denominator 分母多項式の係数
   * @param T 評価する時刻の列
   * @return Gnuplot
   * @throws IOException gnuplotプロセスを起動できない場合
   */
  public static Gnuplot plot(DoubleMatrix numerator, DoubleMatrix denominator, DoubleMatrix T) throws IOException {
    Gnuplot gp = new Gnuplot();
    return plot(gp, numerator, denominator, T);
  }

  /**
   * 単位インパルス応答をプロットする
   * 
   * @param gnuplot gnuplot
   * @param numerator 分子多項式の係数
   * @param denominator 分母多項式の係数
   * @param T 評価する時刻の列
   * @return Gnuplot
   */
  public static Gnuplot plot(Gnuplot gnuplot, DoubleMatrix numerator, DoubleMatrix denominator, DoubleMatrix T) {
    DoubleRationalPolynomial g = Tf2tfn.tf2tfn(numerator, denominator);
    return plot(gnuplot, g, T);
  }

  /**
   * @param G 伝達関数行列
   * @param inputNumber 入力番号
   * @param T 評価する時刻の列
   * @return Gnuplot
   * @throws IOException gnuplotプロセスを起動できない場合
   */
  public static Gnuplot plot(DoubleRationalPolynomialMatrix G, int inputNumber, DoubleMatrix T) throws IOException {
    Gnuplot gp = new Gnuplot();
    return plot(gp, G, inputNumber, T);
  }

  /**
   * 単位インパルス応答をプロットする
   * 
   * @param gnuplot gnuplot
   * @param G 伝達関数行列
   * @param inputNumber 入力番号
   * @param T 評価する時刻の列
   * @return Gnuplot
   */
  @SuppressWarnings("nls")
  public static Gnuplot plot(Gnuplot gnuplot, DoubleRationalPolynomialMatrix G, int inputNumber, DoubleMatrix T) {
    int outputSize = G.getRowSize();
    List<DoubleMatrix> yx = impulse(G, inputNumber, T);
    DoubleMatrix Y = yx.get(0);
    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;
      }
    }

    Canvas canvas = gnuplot.createCanvas();
    canvas.setGridVisible(true);
    canvas.setXLabel("t [sec]");
    canvas.setYLabel("y");
    canvas.plot(T, Y, yy);
    return gnuplot;
  }

}