高等电力网络分析.doc

《高等电力网络分析》

—— IEEE30节点电力网络分析

专 业：电力电子与电力传动
同组成员：杨珊 宋晓英 孙长如
导 师：王艳松 马文忠

二〇一二年四月

第一、二章
第一部分：
本组选用IEEE30节点作为分析对象，首先，根据标准数据，画出电力网络图，如图1所示。然后根据网路图，本单元计算了网络的关联矩阵、节点导纳和节点阻抗矩阵以及添加和移去一条支路的处理。

图1 IEEE30节点电力网络图
一、计算关联矩阵：
为了计算关联矩阵，首先对网络进行节点和支路进行编号和标注方向，尤其是道-支关联矩阵，要求支路必须有方向。选取树枝和连枝，重新编号，如图2所示。

图2 有向图
利用Matlab编程，可直接求出节-支关联矩阵A:

然后根据关联矩阵之间的关系，可分别求出回-支关联矩阵、割-支关联矩阵和道-支关联矩阵。
1. 回-支关联矩阵B:

和A的关系：

2.割-支关联矩阵Q：

和A的关系：

3.道-支关联矩阵T：

和A的关系：

具体程序如下：

function IEEE30
[x,y]=xlsread('C:\Documents and Settings\Administrator\work\30节点数据.xls','sheet3','A2:C51');
A=zeros(30,50);A1=zeros(31,50);
for s=1:50
start=x(s,2);
tail=x(s,3);
zong=x(s,1);
A1(start,zong)=1;
A1(tail,zong)=-1;
end
%去掉参考节点的最后一行，降阶
for s=1:30
for j=1:50
A(s,j)=A1(s,j);
end
end
fprintf(‘节-支关联矩阵A=%8.5f\n')
A
for s=1:30
for j=1:30
AT(s,j)=A(s,j); %树支
end
end
for s=1:30
for j=31:50
AL(s,j-30)=A(s,j); %连支
end
end
BL=eye(20);QT=eye(30);
BT=-1*(AL')*inv(AT');
B=[BT,BL];
fprintf('回-支关联矩阵B=%8.5f\n')
B
QL=-BT';
Q=[QT,QL];
fprintf('割-支关联矩阵=%8.5f\n')
Q
TT=(inv(AT))';
TL=zeros(30,20);
T=[TT,TL];
fprintf('道-支关联矩阵T=%8.5f\n')
T
运行结果如下：
A=

B=

Q=

T=

二、计算节点导纳（阻抗）矩阵
在本节中，本组采用了两种方法对网络进行求解节点导纳矩阵Y，法一先求解网络的不定导纳矩阵，然后去掉参考节点或者地，形成定导纳矩阵；法二，根据导纳矩阵的定义，利用网络直接求出Y。通常阻抗导纳矩阵有三种求解方法，即部分网络法、追加树枝支路法和追加连枝支路法，但是由于Y满秩，所以可以通过直接求逆得出阻抗导纳矩阵，简单快捷。
下面是两种方法的程序：
方法1：
function IEEE30b
[B,y1]=xlsread('C:\Documents and Settings\Administrator\work\30节点数据.xls', 'sheet1','A2:I44');
n1=30;%节点数
n2=41;%支路数
Y=zeros(n1,n1)创建节点导纳矩阵
for j=1:n1
for m=1:n2
if B(m,2)==j&B(m,3)~=0 %支路首段与节点j相连，并且不是补偿电容支路
if B(m,4)==0 %支路无变压器，是线路支路
Y(j,j)=Y(j,j)+1/(B(m,7)+i*B(m,8))+i*B(m,9)/2;
p=B(m,3);
Y(j,p)=-1/(B(m,7)+i*B(m,8));
else %支路有变压器，是变压器支路B(m,4)==1
k=B(m,6)/B(m,5);
zt=B(m,7)+i*B(m,8);
y=1/(k*zt);
y1=(k-1)/(k*zt);
y2=(1-k)/(k*k*zt);
Y(j,j)=Y(j,j)+y+y1;
p=B(m,3);
Y(j,p)=-y;
end
elseif B(m,2)==j&B(m,3)==0
Y(j,j)=Y(j,j)+i*B(m,9);
elseif B(m,3)==j %支路末端与节点j相连，并且不是补偿支路
if B(m,4)==0 %支路是线路支路
Y(j,j)=Y(j,j)+1/(B(m,7)+i*B(m,8))+i*B(m,9)/2;
p=Bm,2);
Y(j,p)=-1/(B(m,7)+i*B(m,8));
else %支路有变压器，是变压器支路B(m,4)==1
k=B(m,6)/B(m,5);
zt=B(m,7)+i*B(m,8);
y=1/(k*zt);
y1=(k-1)/(k*zt);
y2=(1-k)/(k*k*zt);
Y(j,j)=Y(j,j)+y+y2;
p=B(m,2);
Y(j,p)=-y;
end
end
end
end
fprintf('节点导纳矩阵Y=%8.5f\n')
Y
Z=inv(Y);
fprintf('节点阻抗矩阵Z=%8.5f\n')
Z
方法2：
function IEEE30a
%节点导纳矩阵
yb(1,1)=5.2246-j*15.6467;yb(2,2)=1.2437-j*5.096; yb(3,3)=1.7055-j*5.1974;yb(4,4)=8.1954-j*23.5309;yb(5,5)=1.1360-j*4.7725;yb(6,6)=1.6861-j*5.1165;yb(7,7)=6.4131-j*22.3112;
yb(8,8)=2.954-j*7.4493;yb(9,9)=3.5902-j*11.0261;yb(10,10)=6.2893-j*22.0126;yb(11,11)=-j*4.8077;yb(12,12)=-j*1.7986;yb(13,13)=-j*9.0909;yb(14,14)=-j*4.8077;yb(15,15)=-j*3.9063;yb(16,16)=-j*71429;yb(17,17)=1.5266-j*3.1734;yb(18,18)=3.0954-j*6.0973;yb(19,19)=1.9520-j*4.1044;yb(20,20)=2.4910-j*2.2509;
yb(21,21)=1.8678-j*4.3794;yb(22,22)=1.8077-j*3.6914;yb(23,23)=3.0757-j*6.2188;yb(24,24)=5.8824-j*11.7647;yb(25,25)=1.7848-j*3.9854;yb(26,26)=3.9560-j*10.3174;yb(27,27)=5.1019-j*10.9807;yb(28,28)=2.6193-j*5.4008;yb(29,29)=16.7746-j*34.1277;yb(30,30)=1.9683-j*3.9761;yb(31,31)=2.5405-j*3.9544;yb(32,32)=1.4614-j*2.9892;yb(33,33)=1.3099-j*2.2876;
yb(34,34)=1.2183-j*1.8127;yb(35,35)=1.9693-j*3.7602;yb(36,36)=-j*2.5253;yb(37,37)=0.9955-j*1.8810;yb(38,38)=0.6875-j*1.2940;yb(39,39)=0.9121-j*1.7234;yb(40,40)=1.4440-j*4.5408;yb(41,41)=4.3628-j*15.4636; %普通支路数据
yb(42,42)=j*0.0468;yb(43,43)=j*0.0844;yb(44,44)=j*0.0246;yb(45,45)=j*0.0271;yb(46,46)=j*0.0311;yb(47,47)=j*0.0427;yb(48,48)=j*0.0187;yb(49,49)=j*0.0259;yb(50,50)=j*0.0279;
%由支路阻抗数据，求节-支关联矩阵A
A(1,1)=1;A(1,2)=1;A(2,3)=1;A(3,4)=1;A(2,5)=1;A(2,6)=1;A(4,7)=1;A(5,8)=1;A(6,9)=1;A(6,10)=1;A(9,11)=1;A(6,12)=1;A(9,13)=1;A(9,14)=1;A(12,15)=1;A(12,16)=1;A(12,17)=1;A(12,18)=1;A(12,19)=1;A(14,20)=1;A(16,21)=1;
A(15,22)=1;A(18,23)=1;A(19,24)=1;A(10,25)=1;A(10,26)=1;A(10,27)=1;A(10,28)=1;A(21,29)=1;A(15,30)=1;A(22,31)=1;A(23,32)=1;A(24,33)=1;A(25,34)=1;A(25,35)=1;A(28,36)=1;A(27,37)=1;A(27,38)=1;A(29,39)=1;A(8,40)=1;
%输出节点对应的A元素值
A(6,41)=1;A(1,42)=1;A(2,43)=1;A(3,44)=1;A(4,45)=1;A(5,46)=1;A(6,47)=1;A(7,48)=1;A(8,49)=1;A(28,50)=1;
%对地电容支路
A(2,1)=-1;A(3,2)=-1;A(4,3)=-1;A(4,4)=-1;A(5,5)=-1;A(6,6)=-1;A(6,7)=-1;A(7,8)=-1;A(7,9)1;A(8,10)=-1;A(6,11)=-1;A(10,12)=-1;A(11,13)=-1;A(10,14)=-1;A(4,15)=-1;A(13,16)=-1;A(14,17)=-1;A(15,18)=-1;A(16,19)=-1;A(15,20)=-1;
A(17,21)=-1;A(18,22)=-1;A(19,23)=-1;A(20,24)=-1;A(20,25)=-1;A(17,26)=-1;A(21,27)=-1;A(22,28)=-1;A(22,29)=-1;A(23,30)=-1;A(24,31)=-1;A(24,32)=-1;A(25,33)=-1;A(26,34)=1;A(27,35)=-1;
A(27,36)=-1;A(29,37)=-1;A(30,38)=-1;A(30,39)=-1;A(28,40)=-1;A(28,41)=-1;
Y=A*yb*A';
% 6条变压器支路按变比发生变化处理，设原变比为1，变化后变比为t：1，即变压器在原边侧
T1=[9 6 1.0155]; %变压器支路1的首节点、末节点及变比
T2=[6 10 0.9629]; %变压器支路2的首节点、末节点及变比
T3=[12 4 1.0129]; %变压器支路3的首节点、末节点及变比
T4=[28 27 0.9581];% 变压器支路4的首节点末节点及变比
Y(T1(1,1),T1(1,2))=Y(T1(1,1),T1(1,2))+(1-1/T1(1,3))*yb(11,11);Y(T1(1,2),T1(1,1))=Y(T1(1,1),T1(1,2)); % Yji'=Yij'=Yij+detaYij
Y(T1(1,1),T1(1,1))=Y(T1(1,1),T1(1,1))+(1/(T1(1,3)^2)-1)*yb(11,11); % Yii'=Yii+detaYii
Y(T2(1,1),T2(1,2))=Y(T2(1,1),T2(1,2))+1-1/T2(1,3))*yb(12,12);Y(T2(1,2),T2(1,1))=Y(T2(1,1),T2(1,2));
Y(T2(1,1),T2(1,1))=Y(T2(1,1),T2(1,1))+(1/(T2(1,3)^2)-1)*yb(12,12);
Y(T3(1,1),T3(1,2))=Y(T3(1,1),T3(1,2))+(1-1/T3(1,3))*yb(15,15);Y(T3(1,2),T3(1,1))=Y(T3(1,1),T3(1,2));
Y(T3(1,1),T3(1,1))=Y(T3(1,1),T3(1,1))+(1/(T3(1,3)^2)-1)*yb(15,15);
Y(T4(1,1),T4(1,2))=Y(T4(1,1),T4(1,2))+(1-1/T4(1,3))*yb(36,36);Y(T4(1,2),T4(1,1))=Y(T4(1,1),T4(1,2));
Y(T4(1,1),T4(1,1))=Y(T4(1,1),T4(1,1))+(1/(T4(1,3)^2)-1)*yb(36,36);
%2条并联电容支路按支路的添加来修改Y矩阵
l1=[10 0 j*0.19]; %要添加支路1的首末节点及对应的支路导纳
l2=[24 0 j*0.04]; %要添加支路2的首末节点及对应的支路导纳
[m,n]=size(Y);
M1=zeros(m,1); M2=zeros(m,1);
if l1(1,1)~=0
M1(l1(1,1),1)=1;
end
if l1(1,2)~=0
M1(l1(1,2),1)=-1;
end
if l2(1,1)~=0
M2(l2(1,1),1)=1;
end
if l2(1,2)~=0
M2(l2(1,2),1)=-1;
end
Y=Y+M1*l1(1,3)*M1';Y=Y+M2*l2(1,3)*M2'; %Y'=Y-Ml*yl*Ml'
%最后的节点导纳矩阵:
display('节点导纳矩阵如下所示')
Y
Z=inv(Y);
display('节点阻抗矩阵如下所示')
Z
计算结果：

三、添加和移去一条支路
（一）添加一条支路：
相当于在原网络的支路l上并联一个yi的支路

本组选择在原网络上添加支路9-11
（二）移去一条支路：相当于在原网络的支路l上并联一个yi的支路

本组选择在原网络上添加支路6-10
程序如下：
%//////////////////////移去支路6-10////////////////////////////////
[m,n]=size(Y);
l3=[6 10 -5.2246 +i*15.64678];%要移去的支路的首末节点及对应支路导纳
M3=zeros(m,1);
if l3(1,1)~=0
M3(l3(1,1),1)=1;
end
if l3(1,2)~=0
M3(l3(1,2),1)=-1;
end
Y=Y-M3*l3(1,3)*M3';
display('移去支路后新的节点导纳矩阵如下所示：')
Y
计算结果：

%//////////////////////添加支路9-11//////////////////////////////
l4=[9 11 i*9.0909];%要添加支路首末节点及对于支路导纳
M4=zeros(m,1);
if l4(1,1)~=0
M4(l4(1,1),1)=1;
end
if l4(1,2)~=0
M4(l4(1,2),1)=-1;
end
Y=Y+M4*l4(1,3)*M4';
display('添加支路后新的节点导纳矩阵如下所示：')
Y
计算结果：

第三章
本节主要做的内容有：节点优化编号；三角检索存储；Y因子分解的两种方法；求解稀疏线性方程组；Y因子分解图论描述。
一、节点优化编号
1、Tinney-1编号方法：
具体操作：按出线度由小到大进行编号，相同出线度节点先后随意。
优缺点： 简单但编号效果较差。 右图：Tinney-1方法编号Y有向图。
A 、 Tinney_1优化编号实现：
du=zeros(1,30);
for m=1:30 %统计节点出线度
for n=1:30
du(m)=du(m)+Y0(m,n);

end
end
du

min=1;
Min=du(1);
Du=zeros(2,30);
for n=1:30
for m=1:30 %取最小出线度节点
if du(m)
end
end
end
%Du

%------新节点编号--------
Y1=zeros(n1,n1);
for m=1:30
for n=1:30
Y2(m,n)=Y(Du(1,m),Du(1,n))
end
end
B、优化后Y的非零元分布显示：
x=[1:30];y=[1:30];
he1=1;mnx=0;
for m=1:30
for n=1:30
if abs(Y1(m,n))>0
Y0(m,n)=1;
mnx=mnx+1;
x(he1)=m;
y(he1)=n;
he1=he1+1;
end
end
end
mnx
figure (2);
axis([1 30 1 30])
scatter(x,y,'k.')
采用非零元分布显示更加直观了解Y矩阵原始矩阵出线度及节点关系。
C、Y的非零元分布图：

优化前Y非零元分布图

　优化后Y非零元分布图
2、Tinney-2编号方法 (又称半动态节点优化编号方法）
具体操作：按最小出线度编号，编号过程中及时排除已编节点发出边 对未编节点出线度的影响，新生成的边计入后 续节点的出现度。
优缺点：方法较简单，程序复杂度增加不多，效果较好。
本部分编号为手动计算得。
右图：Tinney-2 方法编号法Y图
A、 优化后Y的非零元分布图：

二、Y三角检索存储
1、Y的存储格式：
1、散居格式——修改灵活，检索困难；
2、按行（列）存储格式——修改困难，检索灵活；
3、三角检索存储格式——适合确定的稀疏矩阵；
4、链表存储格式——修改储存灵活。
2、三角检索存储格式：
n U—按行存储Y的上三角部分的非零元素的值
n JU—按行存储Y的上三角部分的非零元素的列号
n IU—存Y中上三角部分每行第一个非零元的序号（行首地址）
n L—按行存储Y的下三角部分的非零元素的值
n IL—按列存储Y的下三角部分的非零元素的行号
n JL—按行存储Y的下三角部分的非零元素的行号
n D—Y的对角元素的值，检索下表不需要存储
3、编程实现Y的三角存储：

[m,n]=size(Y);geshu=0;
for i=1:m
D(1,i)=Y(i,i);
end
IU(1,1)=1;time(1,1)=0;
for i=1:m
cishu=0;
for j=i+1:m
if Y(i,j)~=0
geshu=geshu+1;
cishu=cishu+1;
U(1,geshu)=Y(i,j);
JU(1,geshu)=j;
end
end
time(1,i+1)=cishu;
end
[size2,size1]=size(U);
for i=2:m
IU(1,i)=IU(1,i-1)+time(1,i);
end
IU(1,size1+1)=size1+1;
JL(1,1)=1;time1(1,1)=0;geshu=0;
for j=1:m
cishu1=0;
for i=j+1:m
if Y(i,j)~=0
geshu=geshu+1;
cishu1=cishu1+1;
L(1,geshu)=Y(i,j);
IL(1,geshu)=i;
end
end
time1(1,j+1)=cishu1;
end
[size4,size3]=size(L);
for i=2:1:n
JL(1,i)=JL(1,i-1)+time1(1,i);
end
JL(1,size3+1)=size3+1;

三、Y的LDU分解
方法1：
数值分析的方法—直接按规格化运算后，再进行消去运算，采用全部数据都循环的方法编程。
方法2：
三角检索存储基础上因子分解—在排零存储的基础上，只取非零元U(k)和l(k)进行计算，省去非零元的判断，进行稀疏矩阵的因子分解。
1、因子分解Y->Y’(L*D*U)两种方法对比:
由Y’非零元分布图可以看出，两种方法果是一致的，间接证明程序正确性。
2、优化编号的效果
原始矩阵Y因子分解 Tinney_1优化后Y1因子分解
两图对比易知，节点优化编号大大增大了U矩阵的稀疏度。
3、两种优化编号效果对比

Tinney_1优化后Y1因子分解 Tinney_2优化后Y2因子分解
由程序结果知，图上总出线度为144：140，Tinney_2结果上要优于Tinney_1节点优化编号。
四、利用因子表求解方程组

具体实现：

1、 前代过程
b=I;z=b;
for i=1:n-1
for j=i+1:n
z(j,1)=z(j,1)L(j,i)*z(i,1);

end
end
2、除法运算
for i=1:n
y(j,1)=z(j,1)/D(i,i);
end
3、回代过程
x=y;
for j=n:-1:2
for i=j-1:-1:1
x(i,1)=x(i,1)-U(i,j)*x(j,1);
end
end

五、Y的图论分析
导纳矩阵的稀疏结构可利用图来描述
1、利用图来描述导纳矩阵Y
2、利用图来描述导纳矩阵的因子矩阵U
3、参考习题3.10，分析因子分解过程

右图为Tinney-1优化编号方法下
因子分解后U的赋权有向因子图。

第四章
第四章主要内容：
（一）补偿法求网络方程的修正解
1、前补偿
2、中补偿
3、后补偿
（二）因子表的修正算法
1、因子表的秩1修正（阶次不变）
2、系数矩阵阶次变化时的因子表修正
3、因子表的局部再分解
一、补偿法求网络方程的修正解
1、矩阵求逆辅助定理：
若令nxn阶非奇异矩阵A发生如下变化：
式中，M、N为nxm阶矩阵，a为mxm阶非奇异矩阵，且m≤n,则有：，其条件是可逆。
如果已知A的逆，则可根据上式利用直接求解出的逆。
2、补偿法网络方程的计算：
N维电力系统网络方程为：
当结构或参数发生微小变化而节点注入电流不变时，新的网络方程可写为：，即：，利用矩阵求逆辅助定理有：，其中，（可逆）。
方法一： 后补偿
，或，
方法二：前补偿
或，，。
方法三：中补偿
，，，
三种方法计算量的比较：（m=1）
v 后补偿：FFS+FS+2BS
v 前补偿：FFS+FS+2BS
v 中补偿：2FFS+FS+BS
由上可以看出：中补偿计算量最小。注：补偿法适合注入电流不变，而又要对网络中不同部位发生局部变化时求网络方程的解的应用场合。
3、补偿法求IEEE30节点网络电压方程
，其中。可以根据已知功率S求解，上次我们已经利用系数矩阵的方法求解了该系数网络方程，下面用三种补偿法求节点电压U。
程序：
function zhang3()
load('D:\Program Files\matlab\work\Y.mat');
%///////////////////////LDU分解///////////////////////////////
[m,n]=size(Y);
D=eye(n);
L=eye(n);
U=eye(n);
for q=1:n
for p=1:q-1
U(p,q)=Y(p,q);
for k=1:p-1
U(p,q)=U(p,q)-D(k,k)*L(p,k)*U(k,q);
end
U(p,q)=U(p,q)/D(p,p);
end
D(q,q)=Y(q,q);
for k=1:q-1
D(q,q)=D(q,q)-D(k,k)*L(q,k)*U(k,q);
end
for p=q+1:n
L(p,q)=Y(p,q);
for k=1:q-1
L(p,q)=L(p,q)-D(k,k)*L(p,k)*U(k,q);
end

L(p,q)=L(p,q)/D(q,q);
end
end
%///////////////////////LDU分解///////////////////////////////
[x,x1]=xlsread('D:\shuju\节点注入电流.xls','sheet1','A2:G31');
for m=1:30
U(m,1)=x(m,2)*cos(x(m,3)*3.1415926/180)+i*x(m,2)*sin(x(m,3)*3.1415926/180);
S(m,1)=(x(m,4)+i*x(m,5)-x(m,6)-i*x(m,7))/100;
I1(1,m)=S(m,1)/U(m,1)/sqrt(3);
end
for m=1:30
S1(m,1)=U(m,1)*I1(1,m);
end
I=I1';
% %////////////////////////////前代过程////////////////////
% n=30;
% b=I;z=b;
% for p=1:n-1
% for q=p+1:n
% z(q,1)=z(q,1)-L(q,p)*z(p,1);
% end
% end
% %////////////////////////////除法运算////////////////////
% for q=1:n
% y(q,1)=z(q,1)/D(q,q);
% end
% %////////////////////////////回代运算////////////////////
% x=y;
% for q=n:-1:2
% for m=q-1:-1:1
% x(m,1)=x(m,1)-U(m,q)*x(q,1);
% end
% end
% U=x;
M=zeros(30,1);M(3,1)=1;M(26,1)=-1;yc=10;
% %///////////////////////////网络方程的修正算法—前补偿法
V1=inv(Y)*I;%Ô­À´µÄµçѹ
yy=inv(Y)*M;
c=inv(inv(yc)+M.'*yy)
dV=-1*yy*c*M.'*V1;
V=V1+dV;
% %///////////////////////网络方程的修正算法—中补偿法
yy=inv(Y)*M;
dI=-1*M*c*M.'*inv(Y)*I;
I=dI+I;
V2=inv(Y)*I;
% %///////////////////////////网络方程的修正算法—后补偿法·¨////////////
[xx1]=xlsread('D:\shuju\½Úµã×¢ÈëµçÁ÷.xls','sheet1','A2:G31');
for m=1:30
U(m,1)=x(m,2)*cos(x(m,3)*3.1415926/180)+i*x(m,2)*sin(x(m,3)*3.1415926/180);
S(m,1)=(x(m,4)+i*x(m,5)-x(m,6)-i*x(m,7))/100;
I1(1,m)=S(m,1)/U(m,1)/sqrt(3);
end
for m=1:30
S1(m,1)=U(m,1)*I1(1,m);
end
I=I1';
load('D:\Program Files\matlab\work\Y.mat');
[m,n]=size(Y);L1=eye(m);U1=eye(m);
A=Y;
for k=1:m-1
for p=k+1:m
A(k,p)=A(k,p)/A(k,k);
for q=k+1:n
A(q,p)=A(q,p)-A(q,k)*A(k,p);
end
end
end
for p=1:n
for q=1:p
L1(p,q)=A(p,q);
end
end
for p=1:n
for q=p+1:n
U1(p,q)=A(p,q);
end
end
W=inv(L1)*M;
W1=inv(U1.')*M;
F=inv(L1)*I;
c=inv(inv(yc)+W.'*W1)
dF=-W*c*W1.'*F;
F=F+dF;
V3=inv(U1)*F;
% ///////////////LU·Ö½â///////////////////////////////////////
结果：
LDU分解结果：
矩阵L:

矩阵D:

矩阵U:

前补偿、中补偿、后补偿的结果如下：

二、因子表的修正算法：
1、因子表的秩1修正（阶次不变）
原理：令原来的网络方程系数矩阵Y的因子分解形式为：，网络结构或参数局部发生变化后的系数矩阵变成：，新的矩阵的因子表为：，即：。
程序：%///////////////因子表的秩1修正·¨///////////////////////////////
M=zeros(30,1);M(19,1)=-1;M(21,1)=-1;yc=5;
% a=-inv(yc);
a=yc;
n=30;
for i=1:n-1
A=a*M(i);
D(i,i)=D(i,i)+A*M(i);
B=A/D(i,i);
for j=1:n-i
Mi(j)=M(j+i);
Ui(j)=U(i,j+i);
end
Mi=Mi-Ui*M(i);
Ui=Ui+B*Mi;
% Ui=Ui;
for j=1:n-i
M(j+i)=Mi(j);
U(i,j+i)=Ui(j);
end
A=A-A*B;
end
A=a*M(n);
D(