matlab基于线性准则的考虑风力发电不确定性的分布鲁棒优化机组组合

%% 确定性优化结果 clc close all clear all % A Benchmark Case of Optimal Recourse Under Wind Power Uncertainty D = 320; %MW % demand what = 60; % mean wind output平均风功率 V = 3000; % $/MW % penalty cost for load loss失负荷惩罚 a = 3; % $/MW 发电成本系数 b = 30; % $ % v为预测误差 x = binvar(1); % 机组启停机 p = sdpvar(1); % 机组出力 l = sdpvar(1);% 失负荷 % 不确定性变量 ResultP = []; ResultQ = []; i= 1; for v = -50 :10 :100 Q = b * x + a * p +V * l; constrains = [ p >= 20*x; p <= 300*x ; p + l == D - what - v; l >= 0 ;]; opt = sdpsettings('verbose',1,'solver','cplex'); result = optimize(constrains,Q,opt); % 主问题求解 ResultP(i) = value(p); ResultQ(i) = value(Q); i = i+1; end % % hold on % bar(ResultQ) % ylim([120,320]) xticklabels(-50:10:100) clearvars -except ResultP %% 线性近似求解 %% 模糊集参数 v = sdpvar; u = sdpvar; [A,b]= FuzzyMatrix(); % A*[v;u] <=b; % v,u 无约束 %% 机组相关变量 x = binvar(1); % 一阶段 % 二阶段 p0 = sdpvar; p1 = sdpvar; l0 = sdpvar; l1 = sdpvar; % 相关约束 %% p >= 20*x ( 任意 v,u）1 % min( [p1,0]*[v;u]) >= 20*x - p0 % A[v;u] <=b % v,u 无约束 % 左边转换 DualSize1 = sdpvar(size(A,1),1,'full'); % 对偶变量 DualSize1Constrains = [ A' * DualSize1 == [p1;0]; DualSize1 <= -0.0009 ; b' * DualSize1 >= 20 * x - p0; ]; % 目标函数 %% p <= 300 * x ( 任意 v,u）2 % -p >= -300 * x % min( -p0 - p1 * v - 0 *u ) >=-300 * x % A[v;u] <=b % v,u 无约束 % 左边转换 DualSize2 = sdpvar(size(A,1),1,'full'); % 对偶变量 DualSize2Constrains = [ A' * DualSize2 == [-p1;0]; DualSize2 <= 0; b' * DualSize2 >= -300 * x + p0; ]; % 目标函数 %% p+l = 260 - v ( 任意 v,u）3 % p+l >= 260-v; % p0 + p1*v +0*u + l0 + l1*v + 0*u +v >=260； % min( p1*v +0*u + l1*v + 0*u +v) >= 260 - p0 - l0; % A[v;u] <=b % v,u 无约束 % 左边转换 DualSize3 = sdpvar(size(A,1),1,'full'); % 对偶变量 DualSize3Constrains = [ A' * DualSize3 == [p1+l1+1;0]; DualSize3 <= 0; b' * DualSize3 >= 260-p0-l0; ]; % 目标函数 % p+l <= 260 -v; % -p -l >= v - 260 % -p0 - p1*v - 0*u -l0 - l1*v - 0*u - v >= -260 % min(-p0 - p1*v - 0*u -l0 - l1*v - 0*u - v ) >= -260 % A[v;u] <=b % v,u 无约束 % 左边转换 DualSize4 = sdpvar(size(A,1),1,'full'); % 对偶变量 DualSize4Constrains = [ A' * DualSize4 == [-p1-l1-1;0]; DualSize4 <=0; b' * DualSize4 >= -260 + p0 + l0 ; ]; % 目标函数 %% l >= 0 ( 任意 v,u）4 % l0 + l1*v + 0*u >=0 % min(l0 + l1*v + 0*u) >=0 % A[v;u] <=b % v,u 无约束 % 左边转换 DualSize5 = sdpvar(size(A,1),1,'full'); % 对偶变量 DualSize5Constrains = [ A' * DualSize5 == [l1;0]; DualSize5 <= 0; b' * DualSize5 >= -l0 ; ]; % 目标函数 %% 概率约束相关变量 % dual1 + v * dual2 + u * dual3 >= 3*p + 3000 * l (任意v,u）5 dual1 = sdpvar; dual2 = sdpvar; dual3 = sdpvar; GailvVariableConstrains = [dual3>= 0 ]; % dual1 + dual2 * v + dual3 * u -3*(p0+p1*v+0*u) - 3000*(l0+l1*v + 0*u) >=0; % min(dual1 + dual2 * v + dual3 * u -3*(p0+p1*v+0*u) - 3000*(l0+l1*v +0*u)) >= 0; % A[v;u] <=b % v,u 无约束 % 左边转换 DualSize6 = sdpvar(size(A,1),1,'full'); % 对偶变量 DualSize6Constrains = [ A' * DualSize6 == [dual2-3*p1-3000*l1;dual3]; DualSize6 <= 0; b' * DualSize6 >= -dual1 +3*p0+3000*l0 ; ]; % 目标函数 %% 求解 % 目标函数 % min 30 * x + dual1 + 10 * dual3 func = 30 * x + dual1 + 10* dual3; Constrains = [ A*[v;u] <=b; % 模糊集参数约束 GailvVariableConstrains; % dual3 >=0 DualSize6Constrains; % dual1 + v * dual2 + u * dual3 >= 3*p + 3000 * l (任意v,u）5 DualSize1Constrains; % p >= 20*x ( 任意 v,u）1 DualSize2Constrains; % p <= 300 * x ( 任意 v,u）2 DualSize3Constrains; % p+l = 260 - v ( 任意 v,u）3 DualSize4Constrains; % p+l = 260 - v ( 任意 v,u）3 DualSize5Constrains;]; % l >= 0 ( 任意 v,u）4 opt = sdpsettings('solver','cplex') result = optimize(Constrains,func,opt) %% clc value(p0) value(p1) value(l0) value(l1) k = 1; for i = -50:10:100 Linear(k) = value(p0) + value(p1) * i; % Linear(k) = value(l0) + value(l1) * i; k = k+1; end A=value(A); plot(-50:10:100,Linear) hold on plot(-50:10:100,ResultP) legend('线性近似','确定性优化')