< >另附程序：有解释。</P><>程序一和程序二共同完成问题的第一问。
程序三能够对任意的段预报负荷，根据电力市场规则，爬坡速率等因素，给出各机组的出力分配方案。
程序四是一个优化模型，能够自动的对程序三得出的方案进行阻塞性检验，并自动在阻塞发生时，根据线路上潮流的绝对值百分比尽量小原则和经济原则，给出新的调配方案。
程序五是一个函数文件，与程序四的优化有关。
程序一：
%下面的程序是根据每四组变化的数据，通过一元线形回归，判断单个机组对单个线路的影响，然后，进行假设性检验，如果满足，则置xx相应项为0，否则为1
clear all
clc
a=[ 133.02 73 180 80 125 125 81.1 90;
129.63 73 180 80 125 125 81.1 90;
158.77 73 180 80 125 125 81.1 90;
145.32 73 180 80 125 125 81.1 90;
120 78.596 180 80 125 125 81.1 90;
120 75.45 180 80 125 125 81.1 90;
120 90.487 180 80 125 125 81.1 90;
120 83.848 180 80 125 125 81.1 90;
120 73 231.39 80 125 125 81.1 90;
120 73 198.48 80 125 125 81.1 90;
120 73 212.64 80 125 125 81.1 90;
120 73 190.55 80 125 125 81.1 90;
120 73 180 75.857 125 125 81.1 90;
120 73 180 65.958 125 125 81.1 90;
120 73 180 87.258 125 125 81.1 90;
120 73 180 97.824 125 125 81.1 90;
120 73 180 80 150.71 125 81.1 90;
120 73 180 80 141.58 125 81.1 90;
120 73 180 80 132.37 125 81.1 90;
120 73 180 80 156.93 125 81.1 90;
120 73 180 80 125 138.88 81.1 90;
120 73 180 80 125 131.21 81.1 90;
120 73 180 80 125 141.71 81.1 90;
120 73 180 80 125 149.29 81.1 90;
120 73 180 80 125 125 60.582 90;
120 73 180 80 125 125 70.962 90;
120 73 180 80 125 125 64.854 90;
120 73 180 80 125 125 75.529 90;
120 73 180 80 125 125 81.1 104.84;
120 73 180 80 125 125 81.1 111.22;
120 73 180 80 125 125 81.1 98.092;
120 73 180 80 125 125 81.1 120.44];
b=[165.81 140.13 -145.14 118.63 135.37 160.76
165.51 140.25 -144.92 118.7 135.33 159.98
167.93 138.71 -146.91 117.72 135.41 166.81
166.79 139.45 -145.92 118.13 135.41 163.64
164.94 141.5 -143.84 118.43 136.72 157.22
164.8 141.13 -144.07 118.82 136.02 157.5
165.59 143.03 -143.16 117.24 139.66 156.59
165.21 142.28 -143.49 117.96 137.98 156.96
167.43 140.82 -152.26 129.58 132.04 153.6
165.71 140.82 -147.08 122.85 134.21 156.23
166.45 140.82 -149.33 125.75 133.28 155.09
165.23 140.85 -145.82 121.16 134.75 156.77
164.23 140.73 -144.18 119.12 135.57 157.2
163.04 140.34 -144.03 119.31 135.97 156.31
165.54 141.1 -144.32 118.84 135.06 158.26
166.88 141.4 -144.34 118.67 134.67 159.28
164.07 143.03 -140.97 118.75 133.75 158.83
164.27 142.29 -142.15 118.85 134.27 158.37
164.57 141.44 -143.3 119 134.88 158.01
163.89 143.61 -140.25 118.64 133.28 159.12
166.35 139.29 -144.2 119.1 136.33 157.59
165.54 140.14 -144.19 119.09 135.81 157.67
166.75 138.95 -144.17 119.15 136.55 157.59
167.69 138.07 -144.14 119.19 137.11 157.65
162.21 141.21 -144.13 116.03 135.5 154.26
163.54 141 -144.16 117.56 135.44 155.93
162.7 141.14 -144.21 116.74 135.4 154.88
164.06 140.94 -144.18 118.24 135.4 156.68
164.66 142.27 -147.2 120.21 135.28 157.65
164.7 142.94 -148.45 120.68 135.16 157.63
164.67 141.56 -145.88 119.68 135.29 157.61
164.69 143.84 -150.34 121.34 135.12 157.64];
k=0
xx=zeros(6,8);
for i= 1:6
for j = 1 :8
al=ones(4,2);
al(:,2)=a(4*j-3:4*j,j);%产生一个4行2列的矩阵，这里的第一列为1，第二列为矩阵a的变化的那些数据，每4个为一组。
[bb,bint,r,rint,stats]=regress(b(4*j-3:4*j,i),al);%一元线形回归，这里的bb为回归后的系数
if stats(3)>0.05 %显著性检验，stats是用于检验回归模型的统计量，有三个数值，第一个是R2，其中R是相关系数，第二个是F统计量值，第三个是与统计量F对应的概率P，当P<α时拒绝H0，回归模型成立
xx(i,j)=1;%当回归不成立的时候，让这个位置的数为1，否则为0
end
end
end</P><>程序二
%%%%%%%%%%%%%%下面的是通过矩阵xx，然后去掉xx矩阵为1的相应aa那一列回归矩阵的系数，然后进行线形回归，然后再进行显著性检验，得出一组线路关于机组的线形关系。总共要进行6次。
aa=ones(32,9);%产生一个1矩阵的原因是，在回归的时候，第一行是1
aa(:,2:9)=a;
a1=aa(:,1:8);%产生矩阵 a1，用以进行多元回归，求出机组分别对各个线路的线形关系
a2=a1;a2(:,4:8)=aa(:,5:9);%同上
a3=aa(:,1:7); a3(:,7)=aa(:,9);%同上
a4=a1;a4(:,7:8)=aa(:,8:9);%同上
a5=ones(32,6);a5(:,2:6)=aa(:,3:7);%同上
a6=aa(:,1:7);a6(:,7)=aa(:,8);%同上
jianyan=zeros(1,6);
[b1,bint,r,rint,stats]=regress(b(:,1),a1);%同上，调用回归的函数
if stats(3)>0.05 % 进行多元线形的显著性检验，如果，不显著，或者是非线形的模型，则置jianyan（i）为1，下面的类同。
jianyan(1)=1;
end
[b2,bint,r,rint,stats]=regress(b(:,2),a2);
if stats(3)>0.05
jianyan(2)=1;
end
[b3,bint,r,rint,stats]=regress(b(:,3),a3);
if stats(3)>0.05
jianyan(3)=1;
end
[b4,bint,r,rint,stats]=regress(b(:,4),a4);
if stats(3)>0.05
jianyan(4)=1;
end
[b5,bint,r,rint,stats]=regress(b(:,5),a5);
if stats(3)>0.05
jianyan(5)=1;
end
[b6,bint,r,rint,stats]=regress(b(:,6),a6);
if stats(3)>0.05
jianyan(6)=1;
end
程序三
%%%%%%%%%%%%%%%%%下面的程序，是针对第三问的，即通过电力交易规则，给出清算价最小的一组调配，这个程序可以对任意的段预报进行分配
clear all
clc
demand=982.4 %段预报，这里可以随便更改，也就是说，这样的话，这段程序对任意的段预报，都能给出一个方案来 ，如果取demand=1052。4的话，就是第五问的初步分配方案
a.q=[70 0 50 0 0 30 0 0 0 40 % 段容量
30 0 20 8 15 6 2 0 0 8
110 0 40 0 30 0 20 40 0 40
55 5 10 10 10 10 15 0 0 1
75 5 15 0 15 15 0 10 10 10
95 0 10 20 0 15 10 20 0 10
50 15 5 15 10 10 5 10 3 2
70 0 20 0 20 0 20 10 15 5];
zanshi=a.q;
a.p=[-505 0 124 168 210 252 312 330 363 48 %段价
-560 0 182 203 245 300 320 360 410 495
-610 0 152 189 233 258 308 356 415 500
-500 150 170 200 255 302 325 380 435 800
-590 0 116 146 188 215 250 310 396 510
-607 0 159 173 205 252 305 380 405 520
-500 120 180 251 260 306 315 335 348 548
-800 153 183 233 253 283 303 318 400 800];
biao=[120 73 180 80 125 125 81.1 90]; %当前机组出力方案
papo=[2.2 1 3.2 1.3 1.8 2 1.4 1.8]; %爬坡速度
papo=papo*15;
biaozhi=zeros(6,8);
aa.p=a.p;
aa.q=a.q;
for i=1:8
for j=2:10
aa.q(i,j)=aa.q(i,j)+aa.q(i,j-1);%累加段容量
end
end
liang=0;
zu=zeros(1,8);% zu表示各机组出力分配方案
for i= 1:8
for j=1:10
if aa.q(i,j)>=biao(i)+papo(i)
for k= j+1:10
a.p(i,k)=10000;%由爬坡速度，可以确定，这些段都无法达到，所以，让它的价格达到一极大的值，目的是为了在选数据的时候，选中它;
end
elseif aa.q(i,j)<=biao(i)-papo(i)
a.p(i,j)=-800;% 由爬坡速度，这些段也达不到，，实际上的意思，就是必定选了这些段，所以，让它的价格都取最小，表示，在选数据的时候，一定会选到它
liang=liang+a.q(i,j);% 累加当前所有机组的出力值
zu(i)=zu(i)+a.q(i,j);% 累加各机组的处理值
a.q(i,j)=0;% 累加后，置其为 0，以免重复加如
biaozhi(i,j)=1;
end
end
end
k=0;
for i= 1:8 % 将矩阵转化为一个向量，这样易于排序，同时，保存他们的段价，段容量，行和列
for j= 1:10
k=k+1;
x.p(k)=a.p(i,j);
x.q(k)=a.q(i,j);
x.col(k)=j;
x.row(k)=i;
end
end
x.p;
for i =1:80 % 进行冒泡排序
for j= 79:-1:i
if x.p(j+1)<x.p(j)
z=x.p(j+1);x.p(j+1)=x.p(j);x.p(j)=z;
z=x.q(j+1);x.q(j+1)=x.q(j);x.q(j)=z;
z=x.col(j+1);x.col(j+1)=x.col(j);x.col(j)=z;
z=x.row(j+1);x.row(j+1)=x.row(j);x.row(j)=z;
end
end
end
k=0;
fa=liang;
while fa<=demand % 依照价格，从低到高选取它们的段容量，这样是依照分配要经济的原则
k=k+1;
if aa.q(x.row(k),x.col(k))>biao(x.row(k))+papo(x.row(k)) % 这里，稍微复杂点，表示，如果当前的段容量，大于当前方案与爬坡总和，那么，就只能选该段的一部分
x.q(k)=x.q(k)-(aa.q(x.row(k),x.col(k))-biao(x.row(k))-papo(x.row(k)));
end
fa=fa+x.q(k);
end
for i =1:k
if x.q(i)~=0
biaozhi(x.row(i),x.col(i))=1;% 显示是否选中，若选的话，为 1否则为0
end
end
biaozhi;
x.q(k)=x.q(k)-(fa-demand);% 最后一段，可能只取部分
for i= 1:8
for j= 1:k
if x.row(j)==i
zu(i)=zu(i)+x.q(j);% 累加各机组出力
end
end
end
end
程序四
%%%%%%%下面的这段程序是通过不断的强化约束条件，使在各线路超出限制百分比尽量小的前提下，使得电力公司的出钱（包括付的电费和阻塞费用）尽量小。
%%%%%%%这个模型可以对第4问和第五问求解，只需要修改其中的几个数据和条件就可以。
yuding=zu; % yuding表示预定的方案，上面求出的结果是[150.0000 79.0000 180.0000 99.5000 125.0000 140.0000 95.0000 113.9000]
pc=x.p(k);% pc表示清算价，这里，实际等于选择的最后一个x.p
for i= 1:80
if x.p(i)==10000 %理论上是置其为无穷大，但为了计算，用一个极大的值代替
break
end
end
biao=[120 73 180 80 125 125 81.1 90];
su=[2.2 1 3.2 1.3 1.8 2 1.4 1.8]*15;
x0=[120 73 180 80 125 125 81.1 90];%用constr进行非线形规划时的初始值
v1=biao-su; % x 取值的下限
v3=[153.0000 88.0000 228.0000 99.5000 152.0000 155.0000 102.1000 117.0000] %上限
options(13)=1; %表示 函数fun中，一个g函数是等式约束
m=0;
for j= i-1:-1:k
v2=v3;
v2(x.row(j))=aa.q(x.row(j),x.col(j)-1);% 改变约束条件
[xx,options]=constr('fun',x0,options,v1,v2);%调用非线形规划函数
options(8)
if options(8)<=0 % 说明此时调整后，可以不发生阻塞，如果是取预算为1052。8时，去掉这个项，也就是说，不考虑它的最优植一定得小于0，也就是说，可以发生阻塞
%aa.p(x.row(j),x.col(j)-1)
m=m+1;
mubiao(m)=options(8);% 用向量mubiao保存当前约束条件下进行非先行规划时的最忧值
fangan(m, =xx;%保存当前的分配方案
for ii=1:8 %确定当前方案取得的最后一个段号
Mm=0;
for jj=1:10
Mm=zanshi(ii,jj)+Mm;
if(Mm>xx(ii))
n(ii)=jj;
break
end
end
end
%minprice(m)=0;%存贮此时方案的清算价
minprice(m)=0;
for ii= 1:8 %用minprice 保存当前的清算价
if minprice(m)<aa.p(ii,n(ii))
minprice(m)=aa.p(ii,n(ii));
end
end
zubu(m)=0; %存贮所有可以调配的方案取值时的阻塞补偿。
for ii=1:8 %计算发生阻塞时调配费用，由两部分组成，序内容量和序外容量部分
cha=yuding(ii)-xx(ii);
if cha>0 % 这表示预定比实际上的大，此时，是序内容量，下面的是求出补偿
zubu(m)=zubu(m)+pc*cha*(1-xx(ii)/yuding(ii)); %这是由第二问得出的阻塞费用计算规则的公式
else% 这表示预定比实际上的小，此时，是序外容量，下面的是求出补偿
zubu(m)=zubu(m)+xx(ii)*(aa.p(ii,n(ii))-pc)*(1-yuding(ii)/xx(ii));%这是由第二问得出的阻塞费用计算规则的公式
end
end
% m
% pause
end
end
xi=0;
[zongjia,xi]=min(minprice*982.4/4+zubu/4)
%这里保存了2组结果，第一组是在路超出限制比较小的情况下，电力公司的出力价钱也比较小的情况下的出力方案，第二组是在线路超出限制最小的的情况下的方案，
bilv1=mubiao(xi)%线路超出限制的比率
xuanqufangan1= fangan(xi, %调整后的方案，
qingsuanjia1= minprice(xi)% 清算价
buchangjia1=zubu(xi)/4% 由于输电阻塞而补偿的价钱
zongjia %电力公司付的总价
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%下面的程序是对输电阻塞不可避免时采取的方案
[bilv2,wei]=min(mubiao);
bilv2
xuanqufangan2= fangan(wei,
qingsuanjia2=minprice(wei)
buchangjia2=zubu(wei)/4
zongjia2=qingsuanjia2*982.4/4+buchangjia2</P><P>程序五
function[f,g]=fun(x) % 进行非线形优化时的目标函数和约束条件
y1=110.0164 + 0.0831* x(1) + 0.0488*x(2)+ 0.0532 *x(3)+ 0.1200*x(4)+ -0.0252*x(5)+ 0.1224 *x(6)+ 0.1211*x(7);
y2=131.2189 + -0.0545* x(1) + 0.1279 * x(2) + 0.0333* x(4) + 0.0869 * x(5) + -0.1124* x(6) + -0.0189* x(7) + 0.0987* x(8) ;
y3=-108.3922 + -0.0701* x(1) + 0.0604 * x(2) + -0.1571* x(3) + -0.0101 * x(4) + 0.1238 * x(5) + -0.2021* x(8) ;
y4=78.8167 + -0.0355* x(1) + -0.1046 * x(2) + 0.2044 * x(3) + -0.0212* x(4) + -0.0130* x(5) + 0.1467* x(7) + 0.0753* x(8) ;
y5=131.3861 +0.2449* x(2) + -0.0640* x(3) + -0.0409 * x(4) + -0.0644 * x(5) + 0.0715* x(6) ;
y6=120.7930 + 0.2377* x(1) + -0.0604 * x(2) + -0.0780 * x(3) + 0.0929 * x(4) + 0.0468 * x(5) + 0.1662* x(7) ;
f=max([y1/165-1,y2/150-1,y3/(-160)-1 ,y4/155-1,y5/132-1,y6/162-1]);
g(1)=x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)-982.4;% 这里的982。4是指预报的负荷量，对任意数据都适用。如果取它为1052。4的话，就可以得出第五问的答案
g(2)=y1-186.45;
g(3)=y2-177;
g(4)=-y3-174.4;
g(5)=y4-172.05;
g(6)=y5-151.8;
g(7)=y6-184.68;</P> |
发表于 2004-9-21 11:08:49
楼主