Obsf.java

/*
 * $Id: Obsf.java,v 1.23 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.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.matrix.IntMatrix;
import org.mklab.nfc.matrix.NormType;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.DoubleNumber;
import org.mklab.nfc.scalar.RealNumericalScalar;


/**
 * 可観測部分と負可観測部分を分割するクラスです。
 * 
 * <p>Observable-unobservable decomposition
 * 
 * @author koga
 * @version $Revision: 1.23 $
 * @see org.mklab.tool.control.Ctrf
 * @see org.mklab.tool.control.Obsm
 */
public class Obsf {

  /**
   * 可観測部分と不可観測部分に分解したシステム表現を返します。
   * 
   * <p> もし、<code>rank([[A][C]]) &lt; Rows(A)</code>なら、 相似変換<code>z = T x</code>により
   * 
   * <pre><code> Ab = T*A*T# , Bb = T*B , Cb = C*T#.
   * 
   * [Ano|A12] [Bno] Ab = [---+---], Bb = [---], Cb = [0|Co]. [ 0 |Ao ] [Bo ] </code></pre>
   * 
   * のように変換できます。ただし、 <code>(Co,Ao)</code>は可観測であり、
   * 
   * <pre><code> -1 -1 Co(sI-Ao)Bo = C(sI-A)B. </code></pre>
   * 
   * です。
   * 
   * @param A 元のA行列
   * @param B 元のB行列
   * @param C 元のC行列
   * @return {Ab, Bb, Cb, T, K} (分割されたシステム, 変換行列) decompostion
   */
  public static List<DoubleMatrix> obsf(DoubleMatrix A, DoubleMatrix B, DoubleMatrix C) {
    int n = A.getRowSize();
    DoubleNumber unit = A.getElement(1, 1);

    DoubleNumber tolerance = A.norm(NormType.ONE).multiply(n).multiply(unit.getMachineEpsilon());
    return obsf(A, B, C, tolerance);
  }

  /**
   * 可観測部分空間の次数を決定するために許容誤差<code>tol</code> を用いる。
   * 
   * @param A 元のA行列
   * @param B 元のB行列
   * @param C 元のC行列
   * @param tolerance 許容誤差
   * @return {Ab, Bb, Cb, T, K} (分割されたシステム, 変換行列) decomposition
   */
  public static List<DoubleMatrix> obsf(DoubleMatrix A, DoubleMatrix B, DoubleMatrix C, double tolerance) {
    return obsf(A, B, C, new DoubleNumber(tolerance));
  }

  /**
   * 可観測部分空間の次数を決定するために許容誤差<code>tol</code> を用いる。
   * 
   * @param A 元のA行列
   * @param B 元のB行列
   * @param C 元のC行列
   * @param tolerance 許容誤差
   * @return {Ab, Bb, Cb, T, K} (分割されたシステム, 変換行列) decomposition
   */
  public static List<DoubleMatrix> obsf(DoubleMatrix A, DoubleMatrix B, DoubleMatrix C, DoubleNumber tolerance) {
    String msg;
    if ((msg = Abcdchk.abcdchk(A, B, C)).length() > 0) {
      throw new RuntimeException(msg);
    }

    List<DoubleMatrix> tmp = Ctrf.ctrf(A.conjugateTranspose(), C.conjugateTranspose(), B.conjugateTranspose(), tolerance);
    DoubleMatrix Ac = tmp.get(0);
    DoubleMatrix Bc = tmp.get(1);
    DoubleMatrix Cc = tmp.get(2);
    DoubleMatrix T = tmp.get(3);
    DoubleMatrix K = tmp.get(4);

    DoubleMatrix Ab = Ac.conjugateTranspose();
    DoubleMatrix Bb = Cc.conjugateTranspose();
    DoubleMatrix Cb = Bc.conjugateTranspose();

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

  /**
   * 可観測部分空間の次数を決定するために許容誤差<code>tol</code> を用いる。
   * 
   * @param <RS> スカラーの型
   * @param <RM> 行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @param A 元のA行列
   * @param B 元のB行列
   * @param C 元のC行列
   * @param tolerance 許容誤差
   * @return {Ab, Bb, Cb, T, K} (分割されたシステム, 変換行列) decomposition
   */
  public static <RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> ObsfSolution<RS, RM, CS, CM> obsf(
      RM A, RM B, RM C, RS tolerance) {
    String msg;
    if ((msg = Abcdchk.abcdchk(A, B, C)).length() > 0) {
      throw new RuntimeException(msg);
    }

    CtrfSolution<RS, RM, CS, CM> tmp = Ctrf.ctrf(A.conjugateTranspose(), C.conjugateTranspose(), B.conjugateTranspose(), tolerance);
    RM Ac = tmp.A;
    RM Bc = tmp.B;
    RM Cc = tmp.C;
    RM T = tmp.T;
    IntMatrix K = tmp.K;

    RM Ab = Ac.conjugateTranspose();
    RM Bb = Cc.conjugateTranspose();
    RM Cb = Bc.conjugateTranspose();

    return new ObsfSolution<>(Ab, Bb, Cb, T, K);
    //return new ArrayList<>(Arrays.asList(new DoubleMatrix[] {Ab, Bb, Cb, T, K}));
    //return new MatxList(new Object[] {Ab, Bb, Cb, T, K});
  }

}