Canon.java

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

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

import org.mklab.nfc.eig.DoubleEigenSolution;
import org.mklab.nfc.matrix.DoubleComplexMatrix;
import org.mklab.nfc.matrix.DoubleMatrix;


/**
 * 1入力1出力系をその正準系に変換するクラスです。
 * 
 * <p>Canonical form transformation for SISO Systems
 * 
 * @author koga
 * @version $Revision: 1.24 $
 * @see org.mklab.tool.control.Ctrm
 * @see org.mklab.tool.control.Ctrf
 * @see org.mklab.tool.control.Obsm
 * @see org.mklab.tool.control.Obsf
 */
public class Canon {

  /**
   * システム<code>(A,B,C,D)</code> をモード形式の正準形に変換します。
   * 
   * @param A 元のA行列
   * @param B 元のB行列
   * @param C 元のC行列
   * @param D 元のD行列
   * @return {Ab, Bb, Cb, Db, T} (正準形と変換行列)(canonical-form system)
   */
  public static List<DoubleMatrix> canon(DoubleMatrix A, DoubleMatrix B, DoubleMatrix C, DoubleMatrix D) {
    String type = "modal"; //$NON-NLS-1$
    return canon(A, B, C, D, type);
  }

  /**
   * システム<code>(A,B,C,D)</code> を以下のタイプにしたがって正準形に変換します。
   * 
   * <pre> &quot;modal&quot; : システムの極が対角成分に現れるモード形式 </pre>
   * <pre> &quot;observable&quot; : 特性多項式の係数が最右列に現れる可観測正準形</pre>
   * <pre> &quot;controllable&quot; : 特性多項式の係数が最下行に現れる可制御正準形</pre>
   * <code> T</code> は、変換( <code> x = Tz </code> )のための行列です。
   * 
   * @param A 元のA行列
   * @param B 元のB行列
   * @param C 元のC行列
   * @param D 元のD行列
   * @param type 形式
   * @return {Ab, Bb, Cb, Db, T} (正準形と変換行列)(canonical-form system)
   * 
   */
  public static List<DoubleMatrix> canon(DoubleMatrix A, DoubleMatrix B, DoubleMatrix C, DoubleMatrix D, String type) {
    String message;
    if ((message = Abcdchk.abcdchk(A, B, C, D)).length() > 0) {
      throw new IllegalArgumentException(message);
    }

    int n = A.getRowSize();
    int nt = type.length();

    DoubleMatrix Ab = null, Bb = null, Cb = null, Db = null, T = null;

    if (type.equals("modal".substring(0, Math.min(5, nt)))) { //$NON-NLS-1$
      DoubleEigenSolution tmp = A.eigenDecompose();
      DoubleComplexMatrix val = tmp.getValue();
      DoubleComplexMatrix vec = tmp.getVector();
      int k = 1;
      T = A.createZero(vec.getRowSize(), n);
      while (k <= n) {
        if ((val.getElement(k, k)).getImaginaryPart().isZero() == false) {

          T.setSubMatrix(1, vec.getRowSize(), k, k, (vec.getSubMatrix(1, vec.getRowSize(), k, k)).getRealPart());
          T.setSubMatrix(1, T.getRowSize(), k + 1, k + 1, (vec.getSubMatrix(1, vec.getRowSize(), k, k)).getImaginaryPart());
          k = k + 2;
        } else {
          T.setSubMatrix(1, T.getRowSize(), k, k, vec.getSubMatrix(1, vec.getRowSize(), k, k).getRealPart());
          k = k + 1;
        }
      }
      Ab = T.leftDivide(A).multiply(T);
      Bb = T.leftDivide(B);
      Cb = C.multiply(T);
      Db = D;
    } else if (type.equals("controllable".substring(0, Math.min(12, nt)))) { //$NON-NLS-1$
      if (B.length() > 0) {
        DoubleMatrix V = Ctrm.ctrm(A, B.getSubMatrix(1, B.getRowSize(), 1, 1));
        if (V.rank() < n) {
          throw new RuntimeException(Messages.getString("Canon.1")); //$NON-NLS-1$
        }
        DoubleMatrix Vi = V.inverse();
        T = Obsm.obsm(A, Vi.getSubMatrix(n, n, 1, Vi.getColumnSize())).inverse();
      } else {
        T = A.createZero(0, 0);
      }
      Ab = T.leftDivide(A).multiply(T);
      Bb = T.leftDivide(B);
      Cb = C.multiply(T);
      Db = D;
    } else if (type.equals("observable".substring(0, Math.min(10, nt)))) { //$NON-NLS-1$
      if (C.length() > 0) {
        if (Obsm.obsm(A, C.getSubMatrix(1, 1, 1, C.getColumnSize())).rank() < n) {
          throw new RuntimeException(Messages.getString("Canon.2")); //$NON-NLS-1$
        }
      }
      List<DoubleMatrix> tmp = Canon.canon(A.transpose(), C.transpose(), B.transpose(), D.transpose(), "controllable"); //$NON-NLS-1$
      Ab = tmp.get(0).transpose();
      Bb = tmp.get(2).transpose();
      Cb = tmp.get(1).transpose();
      Db = tmp.get(3).transpose();
      T = tmp.get(4).inverse();
    } else {
      throw new RuntimeException(Messages.getString("Canon.4")); //$NON-NLS-1$
    }

    return new ArrayList<>(Arrays.asList(new DoubleMatrix[] {Ab, Bb, Cb, Db, T}));
  }

}