问题

设计思路

先用传统的基荷、腰荷、峰荷策略解决第一问，然后借助第一问的结果解决第二问。最后后将第一问的策略用代码实现，解决第三问。

一.

1. 划分负荷

明显，无论在任何负载率下，A、B、C发电的经济性是依次变差的，再加上C机组的启停都有额外的费用。所以它们出力的优先级也应当是依次递减。

将日负荷曲线分为基荷、腰荷和峰荷。在10%的旋转备用率条件下，只开A、B机组能够达到的开机容量最多能到达10800MW，无法满足所有负荷需求，所以C机组必须作为备用，在峰荷时出力。
使A、B机组满足基荷、腰荷的负荷要求，A、B、C机组满足峰荷的负荷要求。

日负荷直方图

观察日负荷直方图，有两种划分负荷的方案。

方案甲：5000MW以下为基荷，5000~8000MW为腰荷，8000MW以上为峰荷。
方案乙：5000MW以下为基荷，5000~9000MW为腰荷，9000MW以上为峰荷。

方案甲：
约束条件：
式中A、B分别表示A、B电站机组当日开机的数量。

A、B的开机容量必须大于腰荷，然后提供10%的旋转备用率。
$1000A+600B\geq8800$
A、B的最小技术出力必须小于基荷，自动满足10%的旋转备用率。
$1000A\times50 \% + 600 B\times60\%\le5500$
A、B的数量满足条件
$0 \le A \le 6 \\ 0\le B\le 8$
峰荷时，系统的开机容量需满足旋转备用率的要求，假设负荷率为100%时，C机组全开。
$1000A+600B+350\times4\geq10000\times\left(1+10\%\right)$

解不等式组，无解，下面寻找最接近的整数解。
因为A无论在何种条件下都是最经济的选择，所以先假设：
$A=6$
得：
$\left\{\begin {matrix} A=6\\ B=6 \end{matrix}\right.$
此时，不满足约束条件2：
$1000A\times50\%+600B\times60\%=5160\approx5000$
实际系统中，发电量可以略高于用电量。

方案乙：
约束条件：
式中A、B分别表示A、B电站机组当日开机的数量。

A、B的开机容量必须大于腰荷，然后提供10%的旋转备用率。
$1000A+600B\geq9900$
A、B的最小技术出力必须小于基荷，自动满足10%的旋转备用率。
$1000A\times50\%+600B\times60\%\le5000$
A、B的数量满足条件。
$0\le A\le6\\ 0\le B\le8$
峰荷时，系统的开机容量需满足旋转备用率的要求，假设负荷率为100%时，C机组全开。
$1000A+600B+350×4≥10000×(1+10\%)$
解不等式组（1，2，3），无解。

所以选择方案甲：
$\left\{\begin {matrix} A=6\\ B=6 \end{matrix}\right.$

2. 计算A、B运行成本

下面计算基荷和腰荷的各个工况下A、B每台发电机的出力，已解决在一定条件下成本最小的规划问题。
约束条件：
$P=1000X_{a1}+1000X_{a2}+1000X_{a3}+1000X_{a4}+1000X_{a5}+1000X_{a6}+ 600X_{b1}+600X_{b2}+600X_{b3}+600X_{b4}+600X_{b5} \\ 500\le\ X_{a1},X_{a2},X_{a3},X_{a4},X_{a5},X_{a6}\le1000 \\ {360\le X}_{b1},X_{b2},X_{b3},X_{b4},X_{b5}\le600 \\ minF=\sum_{i=1}^{6}{X_{ai}\left(300-10\frac{X_{ai}}{1000}\right)\times\frac{600}{1000}+\sum_{i=1}^{5}{X_{bi}\left(322-12\frac{X_{bi}}{600}\right)\times\frac{600}{1000}}}$
可以发现，此规划问题的目标函数是二次函数，此问题为二次规划问题，无法用线性规划的单纯形法。所以需采取其他方法。
B机组出力的调节范围为：
${minP}_b=600\times5\times60\%=1800\mathrm{MW} \\ {maxP}_b=600\times5\times100\%=3000\mathrm{MW}$
A机组出力的调节范围为：
$minP_a=1000\times6\times50\%=3000\mathrm{MW}\\ maxP_a=1000\times6\times100\%=6000\mathrm{MW}$
由于B机组的单位成本始终高于A机组，所以尽量减小B机组的出力即可减小整体发电成本。于是，可以将基荷和腰荷时期不同的负荷率对应到A、B不同的出力水平，如下：

负荷/MW	5500	6500	8000
A出力功率/MW	3700	4700	6000
B出力功率/MW	1800	1800	2000

使B的出力始终最小即可。

确定了A、B机组的出力之和，下一步即需确定A、B机组的每台出力，即：在总出力确定的情况下，如何为每台机组分配出力，使得总成本最小。

对于A机组：
如某一小时的A机组运行成本
$F_a=\sum_{i=1}^{6}{X_{ai}\left(300-10\frac{X_{ai}}{1000}\right)\times\frac{600}{1000}} =\sum_{i=1}^{6}\left(-\frac{6}{1000}{X_{ai}}^2+180X_{ai}\right)$
$F_a$ 是 $X_{ai}$ 的开口向下的二次函数，所以 $F_a$ 是 $X_{ai}$ 的上凸函数。
由琴生不等式及其推论可知，应使机组的出力尽量集中于最大出力和最小技术出力处，仅有一台机组处于最大出力和最小技术出力之间时，成本最小。

琴生不等式：
当 $a_1+a_2+a_3+\ldots+a_n=1，且f\left(X\right)$ 为凸函数时：
$a_1f\left(X_1\right)+a_2f\left(X_2\right)+a_3f\left(X_3\right)+\ldots+a_nf\left(X_n\right)\le \ f\left(X_1+X_2+X_3+\ldots+X_n\right)$

B机组同理

所以，可以得出基荷和腰荷时所有发电机的出力水平表如下：

负荷	5500	6500	8000
$X_{a1}$	840	1000	1000
$X_{a2}$	500	1000	1000
$X_{a3}$	500	840	1000
$X_{a4}$	500	500	1000
$X_{a5}$	500	500	1000
$X_{a6}$	500	500	840
$X_{b1}$	360	360	360
$X_{b2}$	360	360	360
$X_{b3}$	360	360	360
$X_{b4}$	360	360	360
$X_{b5}$	360	360	360
$X_{b6}$	360	360	360
总成本	997447.2	1168447	1424947

3. 计算C启停及运行成本

下面考虑峰荷时即负荷率为90%和100%的的情况：

负荷率90%时
为了保证10%的旋转备用率，必须使总开机容量到达 $9900\mathrm{MW}$
$1000A+600B+350C=9900\mathrm{MW}$
所以最少开启2台C机组 $350\mathrm{MW}$ ，总开机容量达到 $10600\mathrm{MW}$ 。
C机组出力调整为最小：
$20\%\times350=70\mathrm{MW}$
则C机组总共出力 $210\mathrm{MW}$ ，A、B需出力：
$9000-2\times70=8860\mathrm{MW}$
负荷率100%时
为了保证10%的旋转备用率，必须使总开机容量到达 $11000\mathrm{MW}$ ，
$1000A+600B+350C=11000\mathrm{MW}$
所以最少开启4台C机组 $350\mathrm{MW}$ ，总开机容量达到 $11000\mathrm{MW}$ 。

目标出力	5500	6500	8000	9000	10000
$X_{a1}$	840	1000	1000	1000	1000
$X_{a2}$	500	1000	1000	1000	1000
$X_{a3}$	500	840	1000	1000	1000
$X_{a4}$	500	500	1000	1000	1000
$X_{a5}$	500	500	1000	1000	1000
$X_{a6}$	500	500	1000	1000	1000
$X_{b1}$	360	360	600	600	600
$X_{b2}$	360	360	600	600	600
$X_{b3}$	360	360	600	600	600
$X_{b4}$	360	360	410	410	600
$X_{b5}$	360	360	360	360	600
$X_{b6}$	360	360	360	360	600
$X_{c1}$	0	0	0	70	190
$X_{c2}$	0	0	0	0	70
$X_{c3}$	0	0	0	0	70
$X_{c4}$	0	0	0	0	70
总出力	5500	6500	8000	9000	10000
开机容量	9360	9360	9360	9950	11000
总成本	997447.2	1168447	1424947	1625101	1900940

4. 结论

得出电站日运行规划

时间（h）	目标出力（MW）	各机组出力（MW）	实际出力（MW）	开机容量（MW）	单位小时出力成本（元）	开停机成本（元）
$X_{a1}$	$X_{a2}$	$X_{a3}$	$X_{a4}$	$X_{a5}$	$X_{a6}$	$X_{b1}$	$X_{b2}$	$X_{b3}$	$X_{b4}$	$X_{b5}$	$X_{b6}$	$X_{c1}$	$X_{c2}$	$X_{c3}$	$X_{c4}$
1	6500	1000	1000	840	500	500	500	360	360	360	360	360	360	0	0	0	0	6500	9360	1168447
2	5000	840	500	500	500	500	500	360	360	360	360	360	360	0	0	0	0	5500	9360	997447.2
3	5000	840	500	500	500	500	500	360	360	360	360	360	360	0	0	0	0	5500	9360	997447.2
4	5000	840	500	500	500	500	500	360	360	360	360	360	360	0	0	0	0	5500	9360	997447.2
5	5000	840	500	500	500	500	500	360	360	360	360	360	360	0	0	0	0	5500	9360	997447.2
6	6500	1000	1000	840	500	500	500	360	360	360	360	360	360	0	0	0	0	6500	9360	1168447
7	6500	1000	1000	840	500	500	500	360	360	360	360	360	360	0	0	0	0	6500	9360	1168447
8	8000	1000	1000	1000	1000	1000	840	360	360	360	360	360	360	0	0	0	0	8000	9360	1424947
9	8000	1000	1000	1000	1000	1000	840	360	360	360	360	360	360	0	0	0	0	8000	9360	1424947
10	9000	1000	1000	1000	1000	1000	1000	600	600	600	410	360	360	70	0	0	0	9000	9950	1625101	36000
11	9000	1000	1000	1000	1000	1000	1000	600	600	600	410	360	360	70	0	0	0	9000	9950	1625101
12	8000	1000	1000	1000	1000	1000	840	360	360	360	360	360	360	0	0	0	0	8000	9360	1424947	36000
13	8000	1000	1000	1000	1000	1000	840	360	360	360	360	360	360	0	0	0	0	8000	9360	1424947
14	8000	1000	1000	1000	1000	1000	840	360	360	360	360	360	360	0	0	0	0	8000	9360	1424947
15	8000	1000	1000	1000	1000	1000	840	360	360	360	360	360	360	0	0	0	0	8000	9360	1424947
16	8000	1000	1000	1000	1000	1000	840	360	360	360	360	360	360	0	0	0	0	8000	9360	1424947
17	9000	1000	1000	1000	1000	1000	1000	600	600	600	410	360	360	70	0	0	0	9000	9950	1625101	36000
18	9000	1000	1000	1000	1000	1000	1000	600	600	600	410	360	360	70	0	0	0	9000	9950	1625101
19	10000	1000	1000	1000	1000	1000	1000	600	600	600	600	600	600	190	70	70	70	10000	11000	1900940	108000
20	10000	1000	1000	1000	1000	1000	1000	600	600	600	600	600	600	190	70	70	70	10000	11000	1900940
21	9000	1000	1000	1000	1000	1000	1000	600	600	600	410	360	360	70	0	0	0	9000	9950	1625101	108000
22	8000	1000	1000	1000	1000	1000	840	360	360	360	360	360	360	0	0	0	0	8000	9360	1424947	3600
23	8000	1000	1000	1000	1000	1000	840	360	360	360	360	360	360	0	0	0	0	8000	9360	1424947
24	6500	1000	1000	840	500	500	500	360	360	360	360	360	360	0	0	0	0	6500	9360	1168447
																			成本(元)：	33415488.17	396000
																			总成本(元)：	33811488.17

二.

1.建立停运容量概率模型

对于确切状态概率，有递推公式
$p(X)=(1-r)p'(X)+rp'(X-C)$
$p(X)$ 为新增一台机组（容量 $C\ \mathrm{MW}$ ， $FOR=r$ ）后，系统停运容量为 $X$ 的确切概率，
对于第一台机组， $p(0)=1-r,p(C)=r$
当 $X<C$ 时， $p'(X-C)=0$
累积状态概率公式
$P(X)=(1-r)P'(X)+rP'(X-C)$
取步长为 $50\mathrm{MW}$ 编程迭代计算停运概率。

2.matlab程序代码

\main.m

%main.m
%华中科技大学 电气与电子工程学院
%强电磁与新技术国家重点实验室
%季天泽
%M201877109

clear
clc
cmat=[ones(6,1)*1000;ones(8,1)*600;ones(4,1)*350]';
rmat=[ones(6,1)*0.1;ones(8,1)*0.08;ones(4,1)*0.12]';

p=[];
p(1)=1;

for i=1:18
    p2=[];%p2代表p(x)
    c=cmat(i);
    r=rmat(i);
    p1=p;%p1代表p'(x)
    p1(x2i(sum(cmat(1:i))))=0;
    for x=0:50:sum(cmat(1:i))%x从0到装机容量
        if x<c
            pxc=0;
        else
            pxc=p1(x2i(x-c));
        end
        px=p1(x2i(x));
        p2=[p2,(1-r)*px+r*pxc];
    end
    p=p2;
end
P=p*(tril(ones(245)));

function  i= x2i(x)
%此函数用于将x值转换为i值
i=x/50+1;
end

三.

筛选出所有的可行的装机组合，然后将其对应的待机组合求出。
将“一”中的策略用代码实现，输入待机组合，输出机组出力表和日运行成本，
同时根据装机组合求其装机成本。
比较各种方案的总成本。

程序代码

%main.m

%华中科技大学 电气与电子工程学院
%强电磁与新技术国家重点实验室
%季天泽
%M201877109

clc
clear
global PowerPlantCombine;%存放A、B、C电站的个数，（1x3）矩阵


LoadRate=[65,55,55,55,55,65,65,80,80,90,90,80,80,80,80,80,90,90,100,100,90,80,80,65]/100;
MaxLoad=12500;
Load=LoadRate*MaxLoad;

 %PowerPlantNumberMat=[6;6;4];
PowerPlantCombineMat=candidatePowerPlantCombine();

dayCostmat=[];


for i = 1:length(PowerPlantCombineMat)%遍历所有可能的电厂组合方式
    PowerPlantCombine=PowerPlantCombineMat(:,i);
    
    dayCost=0;
    Outputmat=[];
    
    for j=1:24%遍历24小时，根据负载求每个电厂的出力组合
        
        load=Load(j);
        
        Output=output(load);%出力组合
        Cost=cost(Output);%出力费用
        
                  
        if j==1
            dayCost=dayCost+Cost;
        elseif j==24
            dayCost=dayCost+Cost+CSwitchCost(Outputmat(1,:),Output);%考虑24点开1点关的C机组的开停机费用
        else
            dayCost=dayCost+Cost+CSwitchCost(Outputmat(end,:),Output);%考虑C机组的开停机费用
        end
        
        Outputmat=[Outputmat;Output];%24小时的出力表
        
    end
    
    dayCostmat=[dayCostmat,dayCost];
end

powerPlantBuildCombineMat=convert(PowerPlantCombineMat);%根据待机容量求相应的装机容量。（待机容量即为当日开机了的所有机组的总容量）


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
buildCostMat=([44067.24,29745.387,12333.3]*10000)*(powerPlantBuildCombineMat-PowerPlantCombineMat);%计算装机成本
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


costmat=dayCostmat*365+buildCostMat;%计算总成本

[min,index]=min(costmat);%找到所有费用的最小值



%还需考虑如果最小值不满足C机组的启停约束，需退而求其次

aim=[costmat;powerPlantBuildCombineMat;PowerPlantCombineMat];


aim(:,index)

%cost.m

%华中科技大学 电气与电子工程学院
%强电磁与新技术国家重点实验室
%季天泽
%M201877109

function cost = cost(output)
%此函数输入每小时的output（出力）矩阵，输出当小时的发电成本

global PowerPlantCombine;

Anumber=PowerPlantCombine(1);
Bnumber=PowerPlantCombine(2);
Cnumber=PowerPlantCombine(3);

Aoutput=output(1:Anumber);
Boutput=output(Anumber+1:Anumber+Bnumber);
Coutput=output(Anumber+Bnumber+1:end);

Acost=sum((300-10.*Aoutput./1000).*600.*Aoutput./1000);
Bcost=sum((322-12.*Boutput./1000).*600.*Boutput./1000);
Ccost=sum((266-16.*Coutput./1000).*1800.*Coutput./1000);
    

cost=Acost+Bcost+Ccost;

end

%convert.m

%华中科技大学 电气与电子工程学院
%强电磁与新技术国家重点实验室
%季天泽
%M201877109

function powerPlantbuildCombineMat = convert(powerPlantonCombineMat)
%通过此函数输入待机组合输出对应的装机组合
powerPlantbuildCombineMat=[];
for i=1:length(powerPlantonCombineMat)
    powerPlantOnCombine=powerPlantonCombineMat(:,i);
    
    Anumber=powerPlantOnCombine(1);
    Bnumber=powerPlantOnCombine(2);
    Cnumber=powerPlantOnCombine(3);
    if Anumber<6
        Anumber=6;
    end
    if Bnumber<8
        Bnumber=8;
    end
    if Cnumber<4
        Cnumber=4;
    end
    powerPlantOnLoad=Anumber*1000+Bnumber*600+Cnumber*350;
    if powerPlantOnLoad<15000
        Cnumber=Cnumber+ceil((15000-powerPlantOnLoad)/350);
    end
    powerPlantbuildCombine=[Anumber;Bnumber;Cnumber];
    powerPlantbuildCombineMat=[powerPlantbuildCombineMat,powerPlantbuildCombine];
end

%condidatePowerPlant.m

%华中科技大学 电气与电子工程学院
%强电磁与新技术国家重点实验室
%季天泽
%M201877109

function A = candidatePowerPlantCombine()
%此函数输出可行的候选待机发电机组合
%   此处显示详细说明

T=[];
J=[];
K=[];
for j=1:8
    for k=1:10
        for t=1:7
            if (1000*j+600*k>=11000)&(500*j+360*k<=6875)&(1000*j+600*k>=13750-350*t)
                J=[J,j];
                K=[K,k];
                T=[T,t];
            end
        end
    end
end
A=[J;K;T];
end

%CSwitchCost.m

%华中科技大学 电气与电子工程学院
%强电磁与新技术国家重点实验室
%季天泽
%M201877109

function switchCost = CSwitchCost(outputFormer,outputNow)
%   输入两个相邻两小时的机组出力状态，通过比较C机组启停状态，计算C机组起停的费用
%   此处显示详细说明
global PowerPlantCombine;

Anumber=PowerPlantCombine(1);
Bnumber=PowerPlantCombine(2);
Cnumber=PowerPlantCombine(3);


CoutputFormer=outputFormer(Anumber+Bnumber+1:end);
CoutputNow=outputNow(Anumber+Bnumber+1:end);

CBoolFormer=logical(CoutputFormer);
CBoolNow=logical(CoutputNow);

Switch=(CBoolFormer-CBoolNow).^2;%是否开/关，若为“1”则开或关
switchNumber=sum(Switch);%总开/关次数
switchCost=switchNumber*10*1800;%费用
end

%output.m

%华中科技大学 电气与电子工程学院
%强电磁与新技术国家重点实验室
%季天泽
%M201877109

function output = output(load)
%输入负荷，输出所有电站的出力水平矩阵，也就是“一”中的部分
%将出力大小分为四个范围，以0，M1，M2，M3，M4为界，分别考虑出力情况
%M1~M2为基荷
%M2~M3为腰荷
%M3~M4为峰荷
%其他范围无解

global PowerPlantCombine
Anumber=PowerPlantCombine(1);%A机组的数量
Bnumber=PowerPlantCombine(2);%B机组的数量
Cnumber=PowerPlantCombine(3);%C机组的数量

MaxA=Anumber*1000;%A、B、C所有机组的出力范围
MaxB=Bnumber*600;

MaxC=Cnumber*350;
MinA=Anumber*500;

MinB=Bnumber*360;
MinC=Cnumber*70;

M1=MinA+MinB;%0~M1无解
M2=MaxA+MinB;%M1~M2,C停机，B最小出力
M3=MaxA+MaxB;%M2~M3，C停机，A最大出力
M4=MaxA+MaxB+MaxC;%M3~M4，C开机，A、B最大出力（B可能有一台机组需减小出力）
%M4以上，无解


%通过总负荷来得到A、Bnumber、C机组分别的总负荷，对于相同机组，使更多的机组处于出力区间的端点处，一台机组处于出力区间之内（琴生不等式）
if (M1<=load)&&(load<M2)
    Aload=load-MinB;
    Bload=MinB;
    AmaxNumber=floor((Aload-500*Anumber)/500);
    AchangableLoad=(Aload-500*Anumber)-500*AmaxNumber+500;
    outputA=[ones(1,AmaxNumber)*1000,[AchangableLoad],ones(1,Anumber-AmaxNumber-1)*500];
    outputB=ones(1,Bnumber)*360;
    outputC=ones(1,Cnumber)*0;
    
elseif (M2<=load)&&(load<M3)
    Aload=MaxA;
    Bload=load-MaxA;
    
    BmaxNumber=floor((Bload-360*Bnumber)/240);
    BchangableLoad=(Bload-360*Bnumber)-240*BmaxNumber+360;
    outputB=[ones(1,BmaxNumber)*600,BchangableLoad,ones(1,Bnumber-BmaxNumber-1)*360];
    outputA=ones(1,Anumber)*1000;
    outputC=ones(1,Cnumber)*0;
    
elseif (M3<=load)&&(load<M4)
    Aload=MaxA;
    Bload=MaxB;
    Cload=load-MaxA-MaxB;
    
    
    if Cload<70%如果A、B都为最大出力将导致C的目标出力小于其可能的最小出力，那么需将C设置为最小出力(一台机为70，其余为0)，重新设置B的出力。
        
        Cload=70;
        Bload=load-Aload-Cload;
        
        BmaxNumber=floor((Bload-360*Bnumber)/240);
        BchangableLoad=(Bload-360*Bnumber)-240*BmaxNumber+360;
        outputB=[ones(1,BmaxNumber)*600,BchangableLoad,ones(1,Bnumber-BmaxNumber-1)*360];
        outputA=ones(1,Anumber)*1000;
        outputC=ones(1,Cnumber)*0;
        outputC(1)=70;
        
%         outputB=ones(1,Bnumber);
%         outputB(end)=Bload-(Bnumber-1)*600;

    else
        
        ConNumber=ceil(Cload/350);%确认C开机的数量
        CmaxNumber=ceil((Cload-70*ConNumber)/280-1);
        CchangableLoad=(Cload-70*ConNumber)-280*CmaxNumber+70;
        outputC=[ones(1,CmaxNumber)*350,CchangableLoad,ones(1,ConNumber-CmaxNumber-1)*70,ones(1,Cnumber-ConNumber)*0];
        outputB=ones(1,Bnumber)*600;
        outputA=ones(1,Anumber)*1000;
    end
else
    '此发电机组合错误'
    PowerPlantCombine
end

output=[outputA,outputB,outputC];


end

你可以随意更改我的代码的任何部分，如果我的思路和代码对你有所帮助，请点赞或评论，这对我非常重要，谢谢

利用matlab解决电力系统规划问题