风蓄联合优化程序

clearclc
%求负荷预测曲线
L1=[7800 7750 7760 7650 7660   7598 7539 7438 ...
    7450 7337 7378 7337 7282  7239 7169  7115 7080 7050 7000 7090 ];
t=1:96;
L2=[7040 7110 7200 7337 7599 7630 7900 8011 8143 8500 8420 8400 8496 8700 ...
    8900  9280 9270 9529 9685 9781 9900 10200 10100 10230 10360 10650  10500];
L3=[10487 9979 9598];
L4=[ 9591 9650 9759 9790 9870 9980 10000 10111 10130 10210 10300  10280 ...
     10401 10490 10470 10520 10600 10720 10670 10800];
 L5=[10770 10660 10500 10420 10287 10189 10000 10076 10100 10070 9810 9650 ...
     9470 9300 9230 9100 8950 8710 8790 8640 8587 8487 8340 8202 8110 8000];
Lsi=[L1 L2 L3 L4 L5]‘;

%求风功率预测平均值曲线
Pvdi=[5057 5100 5200 5220  5294 5320 5370 5400 5432 5300 5050 4899 4841 4899 4940 4870 ...
    5060 4500 4000 3600 3253 3175 3000 2870 2724 2600 2400 2300 2222 2000 1759 1479 ...
    1251 1480 1700 1900 2115 2450 2789 3000  3377 3510 3657 3800 3923 3423 3000 2600 ...
    2208 2700 2500 4000 4376  4300 4355 4400 4470 4520 4800 5080 5130 5190 5247 5286 5299 5310 5321 ...
    5330 5150 5000 4820 4661 4200 3800 3461 3257 3200 3165 3100 3046 3100 3151 3210 ...
    3244 3100 3000 2800 2752 3000 3333 3628 3949 4000 4100 4280 4460 ]‘;
%求风功率预测值区域
Pvi1=0.95*Pvdi;
Pvi2=1.05*Pvdi;
plot（tLsi/1000‘-xr‘tPvdi/1000‘-.k‘tPvi1/1000‘-k‘tPvi2/1000‘-k‘）;
axis（[0100014]）;
xlabel（‘时间间隔/15min‘）;
ylabel（‘功率/GW‘）;
legend（‘负荷预测‘‘风功率预测平均值‘‘风功率预测值上、下限值‘‘Location‘‘Northwest‘）;
%目标函数f1
c1=zeros（961）;
for k1=1:32
    c1（k1）=0.5;
end
for k2=33:88
    c1（k2）=1;
end
for k3=89:96
    c1（k3）=0.5;
end
c2=（-1）*c1;
c3=（-3000）+zeros（1921）;
C=[c1;c1;c2;c3;zeros（1921）]; %目标函数f1的表达式系数矩阵
%线性等式约束系数矩阵Aeq
Aeq1=zeros（1672）;
for k=97:192
    k1=k+96;
    Aeq1（k）=-0.287;
    Aeq1（k1）=0.217;  %针对式（2345）的1个线性等式约束
end
Aeq2=zeros（95672）;
for k=1:95
    aeq2=zeros（1672）;   %单行矩阵放在循环内，不具有记忆作用
    k1=k+480;
    k2=k+481;
    k3=k+288;
    k4=k+384;
    aeq2（k1）=1;
    aeq2（k2）=-1;
    aeq2（k3）=1;
    aeq2（k4）=-1;
    Aeq2（k:）=aeq2;   %针对式（7）的95个线性等式约束
end
Aeq=[Aeq1;Aeq2];
% disp（Aeq）;
%线性等式约束常数矩阵beq
beq=zeros（961）;
%求线性不等式约束系数矩阵A
a1=zeros（1672）;    %单行矩阵放在循环外，具有记忆作用
A1=zeros（95672）;
for k=1:95
    k1=k+96;
    k2=k+192;
    a1（k1）=（-0.287）;
    a1（k2）=0.217;
    A1（k:）=a1;
end                      %针对式（2345）右边不等式
A2=（-1）.*A1;       %针对式（2345）右边不等式
A3=zeros（1672）;
for k=289:480
    A3（k）=1;
end                     %针对式（8）
A4=zeros（96672）;
for k=1:96
    a4=zeros（1672）;   %单行矩阵放在循环内，不具有记忆作用
    k1=k+192;
    k2=k+480;
    a4（k1）=（-1）;
    a4（k2）=200;
    A4（k: