基于Cplex的发电机组经济运行优化(MATLAB实现)

场景设定：
1）假设有N = 3个不同类型的发电厂。 
2）这些工厂中的每一个都具有最大功率容量，并且在开启时具有最小容量。
3）调度时间尺度为24h，时间步长为1h。预测的电力需求由周期函数/随机分布/正态分布给出。
4）一个时间单位运行工厂i的成本由下式二次函数给出:
                                                              Q_i*(P_i)^2 + C_i*P_i     
其中P_i中是来自工厂的输送功率P_i。

代码：

clear;
clc;
%% 初始化参数N  T  P
Nunits = 3;
Horizon = 24;

Pmax = [120;100;80];
Pmin = [10;20;20];

%Pforecast = 150 + 50*sin((1:Horizon)*2*pi/24);%正弦分布的周期函数（预测1）
%Pforecast = 150 + 50*rand(1，Horizon);%随机分布（预测2）
%Pforecast = 160 + 50*normrnd(1,0.2,[1 Horizon]);%负荷正态分布波动（预测3）
Load_space = randn(Horizon,1) * 0.3 + 1;
Load_space = Load_space(Load_space > 0 & Load_space < 2);
Pforecast = 150 + 50*Load_space(1:Horizon);%限定范围的正态分布波动（预测4）

%% 变量格式
%onoff是二进制变量，其控制单元i是否在时间k开启；
%第二个变量是表示在时间k由电厂i输送的功率的连续变量。
onoff = binvar(Nunits,Horizon,'full');
P     = sdpvar(Nunits,Horizon,'full');

%% 约束条件
Constraints = [];
%约束1：任意时刻k，机组出力应当在（min,max)之间
for k = 1:Horizon
    Q = diag([random('Poisson',0.07) 0.08 0.06]);
           %随便设置了其中一个Q为泊松分布，可以随便调整
    C = [20 10 20];
  Constraints = [Constraints, onoff(:,k).*Pmin <= P(:,k) <= onoff(:,k).*Pmax];
                
end
%约束2：任意时刻k，N个机组出力之和应当满足预测的电力需求
for k = 1:Horizon
  Constraints = [Constraints, sum(P(:,k)) >= Pforecast(k)];
end

%% 优化目标
Objective = 0;
for k = 1:Horizon
  Objective = Objective + P(:,k)'*Q*P(:,k) + C*P(:,k);
  %计算预测范围内的总运行成本
end

%% yalmip调用cplex求解器
ops=sdpsettings('solver', 'cplex');
optimize(Constraints,Objective,ops)

%% 绘制阶梯图
stairs(value(P)');
legend('Unit 1','Unit 2','Unit 3');