delphi - delphi中复数(NxN)矩阵的逆矩阵

Question

一般来说，我想计算复方 (NxN) 矩阵的逆矩阵。

F.ex 我有一个 5x5 矩阵:

Ybus = [ 
   6.2500 -18.6950i, -5.0000 +15.0000i, -1.2500 + 3.7500i,   0,                 0 ; 
  -5.0000 +15.0000i, 10.8333 -32.4150i, -1.6667 + 5.0000i,  -1.6667 + 5.0000i, -2.5000 + 7.5000i; 
  -1.2500 + 3.7500i, -1.6667 + 5.0000i, 12.9167 -38.6950i, -10.0000 +30.0000i, 0; 
   0,                -1.6667 + 5.0000i, -10.0000+30.0000i,  12.9167 -38.6950i, -1.2500 + 3.7500i; 
   0,                -2.5000 + 7.5000i, 0,                  -1.2500 + 3.7500i,  3.7500 -11.2100i;
]

如何使用 Delphi 计算该矩阵的逆矩阵 (Zbus = inverse (Ybus))？

Answer 1

存在 Matrix libray for Delphi由 Nikolai Shokhirev 先生(GNU2 许可)编写。它不是一个完整的库，而是一个很好的起点。除此之外，它仅能够计算实值矩阵的逆矩阵。但是有一种方法可以使用实值矩阵来计算复矩阵的逆:

根据this Matlab ressource ，给定复方矩阵 M = A + iB，其逆也是复方矩阵 Z = X + iY，其中 A、B 和 X、Y 均为实矩阵。发现 M^-1 = Z 或(A + iB)^-1 = (A + BA^-1*B)^-1 - i*(B + A*B^-1*A)^-1前提是涉及求逆的矩阵必须是非奇异的。

以下代码使用矩阵库和Matlab引用来求复数Ybus矩阵的逆，也可以一般用于求复数NxN矩阵的逆:

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, DateUtils,

  // from matrix library
  // http://www.shokhirev.com/nikolai/programs/tools/PasMatLib/download.html
  uDynObjAlg,
  uDynArrays,
  uMatTypes;

type
  TForm1 = class(TForm)
    Memo1: TMemo;
    procedure FormCreate(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

// extra matrix utils:
function CMat2Str(AMatrix: ICArr2D): String;
var
  row : Integer;
  col : Integer;
begin
  Result := '';
  for row := 1 to AMatrix.Dim1 do
  begin
    for col := 1 to AMatrix.Dim2 do
    begin
       Result := Result + CmplxToStr0(AMatrix[row, col], 10, 3) + ' ';
    end;
    Result := Result +#13#10;
  end;
  Result := Result;
end;

function Mat2Str(AMatrix: IFArr2D): String;
var
  row : Integer;
  col : Integer;
  s : string;
begin
  Result := '';
  for row := 1 to AMatrix.Dim1 do
  begin
    for col := 1 to AMatrix.Dim2 do
    begin
       Str(AMatrix[row, col]:10:3,s);
       Result := Result + s + ' ';
    end;
    Result := Result + ';'+#13#10;
  end;
end;

function MtAddMt(const M1: IFArr2D; const M2: IFArr2D): IFArr2D;
var
  row: TInt;
  col : TInt;
  t: IFArr2D;
begin
  if (M1.Lo1<>M2.Lo1) or ( M1.Hi1<>M1.Hi1) or (M1.Lo2<>M2.Lo2) or ( M1.Hi2<>M1.Hi2) then
    Raise ERangeError.Create(RS_LimMismatch);
  t := TFArr2D.Create(M1,true);
  for row := t.Lo1 to t.Hi1 do
    for col := t.Lo2 to t.Hi2 do
      t[row,col] := t[row,col] + M2[row, col];
  result := t;
end;


procedure TForm1.FormCreate(Sender: TObject);
const
  cInversionCount = 1000;

var
  YBus : ICArr2D;
  ZBus : ICArr2D;

  YBusRe: IFArr2D;
  YBusIm: IFArr2D;
  YBusReInv : IFArr2D;
  YBusImInv : IFArr2D;

  row, col : Integer;

  ZBusRe : IFArr2D;
  ZBusIm : IFArr2D;

  timeStart: TDateTime;
  timeStop : TDateTime;

  n : Integer;

begin
  YBus := TCArr2D.Create(1,5, 1,5);

  // fill matrix:

  // row 1:
  // 6.2500 -18.6950i, -5.0000 +15.0000i, -1.2500 + 3.7500i, 0, 0 ;
  YBus.Value[1,1] := cmplx(  6.2500, -18.6950 );
  YBus.Value[1,2] := cmplx( -5.0000,  15.0000 );
  YBus.Value[1,3] := cmplx( -1.2500, 3.7500 );
  YBus.Value[1,4] := cmplx( 0, 0 );
  YBus.Value[1,5] := cmplx( 0, 0 );

  // row 2:
  // -5.0000 +15.0000i, 10.8333 -32.4150i, -1.6667 + 5.0000i, -1.6667 + 5.0000i, -2.5000 + 7.5000i;
  YBus.Value[2,1] := cmplx( -5.0000, 15.0000 );
  YBus.Value[2,2] := cmplx( 10.8333,  -32.4150 );
  YBus.Value[2,3] := cmplx( -1.6667, 5.0000 );
  YBus.Value[2,4] := cmplx( -1.6667, 5.0000 );
  YBus.Value[2,5] := cmplx( -2.5000, 7.5000 );

  // row 3:
  // -1.2500 + 3.7500i, -1.6667 + 5.0000i, 12.9167 -38.6950i, -10.0000 +30.0000i, 0;
  YBus.Value[3,1] := cmplx( -1.2500, 3.7500 );
  YBus.Value[3,2] := cmplx( -1.6667,   5.0000 );
  YBus.Value[3,3] := cmplx( 12.9167, -38.6950 );
  YBus.Value[3,4] := cmplx( -10.0000, 30.0000 );
  YBus.Value[3,5] := cmplx( 0, 0 );

  // row 4:
  // 0, -1.6667 + 5.0000i, -10.0000 +30.0000i, 12.9167 -38.6950i, -1.2500 + 3.7500i;
  YBus.Value[4,1] := cmplx( 0, 0 );
  YBus.Value[4,2] := cmplx( -1.6667,   5.0000 );
  YBus.Value[4,3] := cmplx( -10.0000, 30.0000 );
  YBus.Value[4,4] := cmplx( 12.9167, -38.6950 );
  YBus.Value[4,5] := cmplx( -1.2500, 3.7500 );

  // row 5:
  // 0, -2.5000 + 7.5000i, 0, -1.2500 + 3.7500i, 3.7500 -11.2100i
  YBus.Value[5,1] := cmplx( 0, 0 );
  YBus.Value[5,2] := cmplx( -2.5000,   7.5000 );
  YBus.Value[5,3] := cmplx( 0, 0 );
  YBus.Value[5,4] := cmplx( -1.2500, 3.7500 );
  YBus.Value[5,5] := cmplx( 3.7500, -11.2100 );


  // compute inverse of complex matrix using relation:
  // http://www.mathworks.com/matlabcentral/fileexchange/49373-complex-matrix-inversion-by-real-matrix-inversion
  // Given a complex square matrix M = A + i*B,
  // its inverse is also a complex square matrix Z = X + i*Y,
  // where A, B and X, Y are all real matrices. It is found that
  //         M^-1 = Z or
  // (A + i*B)^-1 = (A + B*A^-1*B)^-1 - i*(B + A*B^-1*A)^-1
  // Provided that those matrices involved inversion must be nonsingular.

  // with performance profiling:
  timeStart := now;
  for n := 1 to cInversionCount do
  begin
    // Create real part matrix:
    YBusRe := TFArr2D.Create( YBus.Lo1, YBus.Hi1, YBus.Lo2, YBus.Hi2);
    for row := 1 to YBus.Dim1 do
    begin
      for col := 1 to YBus.Dim2 do
      begin
         YBusRe[row, col] := YBus[row, col].Re;
      end;
    end;
    // Create imaginary part matrix:
    YBusIm := TFArr2D.Create( YBus.Lo1, YBus.Hi1, YBus.Lo2, YBus.Hi2);
    for row := 1 to YBus.Dim1 do
    begin
      for col := 1 to YBus.Dim2 do
      begin
         YBusIm[row, col] := YBus[row, col].Im;
      end;
    end;

    // compute inverse of real matrices:
    YBusReInv := PseudoinverseMt( YBusRe );
    YBusImInv := PseudoinverseMt( YBusIm );

    // compute:
    // (A + B*A^-1*B)^-1 - i*(B + A*B^-1*A)^-1
    ZBusRe := PseudoinverseMt( MtAddMt( YBusRe, MtxMt( MtxMt(YBusIm, YBusReInv ), YBusIm ) ) );
    ZBusIm := PseudoinverseMt( MtAddMt( YBusIm, MtxMt( MtxMt( YBusRe, YBusImInv ), YBusRe ) ) );

    // and finally combine to inverse complex matrix:
    ZBus :=  TCArr2D.Create( YBus, False );
    for row := 1 to ZBus.Dim1 do
    begin
      for col := 1 to ZBus.Dim2 do
      begin
         ZBus[row, col] := cmplx( ZBusRe[row, col], -ZBusIm[row, col] );
      end;
    end;
  end;
  timeStop := now;

  // print results:
  Memo1.Text := 'YBus = ' + #13#10 + CMat2Str( YBus ) + #13#10+
                'YBusRe = ' + #13#10 + Mat2Str( YBusRe ) + #13#10 +
                'YBusReInv = ' + #13#10 + Mat2Str( YBusReInv ) + #13#10 +
                'Verify inverse, I = YBusRe x YBusReInv =' + #13#10 + Mat2Str( MtxMt(YBusRe, YBusReInv ) ) + #13#10 +

                'YBusIm = ' + #13#10 + Mat2Str( YBusIm ) + #13#10 +
                'YBusImInv = ' + #13#10 + Mat2Str( YBusImInv ) + #13#10 +
                'Verify inverse, I = YBusIm x YBusImInv =' + #13#10 + Mat2Str( MtxMt(YBusIm, YBusImInv ) ) + #13#10 +

                'ZBus = ' + #13#10+ CMat2Str( ZBus ) + #13#10+
                'Verify ZBus, I = YBus x ZBus = ' + #13#10+  CMat2Str( CMtxCMt( YBus, ZBus ) ) + #13#10 +
                'Performance: ' + FormatFloat('0.00', MilliSecondsBetween(timeStop, timeStart ) / cInversionCount) + ' ms for 1 inversion. Or ' +
                IntToStr( Round( 1000 / (MilliSecondsBetween(timeStop, timeStart)/cInversionCount))) + ' inversions per second. (Intel i7-4790 CPU @ 3.60GHz)';

end;

end.