delphi 中复 (NxN) 矩阵的逆
Inverse of complex (NxN) matrix in delphi
一般来说，我想计算复方 (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))？
存在 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.
一般来说，我想计算复方 (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))？
存在 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.