1 主要内容
程序主要复现的是《多源动态最优潮流的分布鲁棒优化方法》，针对大规模清洁能源接入电网引起的系统鲁棒性和经济性协调问题，提出含风–光–水–火多种能源的分布鲁棒动态最优潮流模型。采用分布鲁棒优化方法将风光不确定性描述为包含概率分布信息的模糊不确定集。将模糊不确定集构造为一个以风光预测误差经验分布为中心，以 Wasserstein距离为半径的 Wasserstein 球。在满足风光预测误差服从模糊不确定集中极端概率分布情况下最小化运行费用。由于梯级水电厂模型为混合整数模型，为了提高计算效率，将交流潮流近似为解耦线性潮流,以118节点系统为例验证了方法的可行性。
在研究过程中，有些地方描述不是很清楚，然后参考同作者英文文献获得启发《Wasserstein Metric Based Distributionally Robust Approximate Framework for Unit Commitment》,两个文章公式部分基本一致，表达方式有所差别，所以单一个文献看不明白的话可以找到作者其他文献就更容易了解作者的意图，尤其是有些期刊论文说的比较简约，但是在硕博论文中就解释的非常仔细，一点小心得分享给大家。
2 代码
%% 决策变量 x_theta = sdpvar(nbus, Horizon,'full');%网络角度 V = sdpvar(nbus, Horizon,'full');%网络节点电压 x_P_h = sdpvar(ngen, Horizon,'full');%风光调整前火电 x_P_s = sdpvar(ns, Horizon,'full');%风光调整前水电 % x_P_hz = sdpvar(ngen, Horizon,'full'); % x_P_sz = sdpvar(ns, Horizon,'full'); x_P_w = sdpvar(nw, Horizon,'full'); x_P_v = sdpvar(nv, Horizon,'full'); ww = sdpvar(1,Horizon,'full');%风力偏差 wp = sdpvar(1,Horizon,'full');%光伏偏差 wwp = sdpvar(1,Horizon,'full');%风光总偏差 alfah = sdpvar(ngen,Horizon,'full');%火电机组参与因子 alfas = sdpvar(ns,Horizon,'full');%水电机组参与因子 rgmax = 50.*ones(ngen, Horizon);%火电旋转备用容量 rgmin = 10.*ones(ngen, Horizon);%火电旋转备用容量 rsmax = 50.*ones(ns, Horizon);%水电旋转备用容量 rsmin = 10.*ones(ns, Horizon);%水电旋转备用容量 rhog = 2.*ones(ngen, Horizon); rhos = 3.*ones(ns, Horizon); k1 = sdpvar(1);%对偶变量 k2 = sdpvar(1);%对偶变量 tk1 = sdpvar(1,K,'full');%辅助变量 tk2 = sdpvar(1,K,'full');%辅助变量 %平方分段线性化参数 gn=5;%分段数 x_pf=sdpvar(ngen, Horizon,'full');%p的平方 x=sdpvar(ngen, Horizon,'full'); gw1=sdpvar(gn+1,Horizon,'full');%辅助参数，下同 gw2=sdpvar(gn+1,Horizon,'full'); gw3=sdpvar(gn+1,Horizon,'full'); gw4=sdpvar(gn+1,Horizon,'full'); gw5=sdpvar(gn+1,Horizon,'full'); gw6=sdpvar(gn+1,Horizon,'full'); gw7=sdpvar(gn+1,Horizon,'full'); gw8=sdpvar(gn+1,Horizon,'full'); gw9=sdpvar(gn+1,Horizon,'full'); gw10=sdpvar(gn+1,Horizon,'full'); gw11=sdpvar(gn+1,Horizon,'full'); gw12=sdpvar(gn+1,Horizon,'full'); gw13=sdpvar(gn+1,Horizon,'full'); gw14=sdpvar(gn+1,Horizon,'full'); gz1=binvar(gn, Horizon,'full');gz2=binvar(gn, Horizon,'full');gz3=binvar(gn, Horizon,'full');gz4=binvar(gn, Horizon,'full');gz5=binvar(gn, Horizon,'full'); gz6=binvar(gn, Horizon,'full');gz7=binvar(gn, Horizon,'full');gz8=binvar(gn, Horizon,'full');gz9=binvar(gn, Horizon,'full');gz10=binvar(gn, Horizon,'full'); gz11=binvar(gn, Horizon,'full');gz12=binvar(gn, Horizon,'full');gz13=binvar(gn, Horizon,'full');gz14=binvar(gn, Horizon,'full'); %% 约束条件生成 cons = []; cons = [cons,wwp == ww + wp]; % 火电 Phmax = 10.*[460;300;443;320;330;460;300;443;320;330;460;300;443;320];%火电机组上限 Phmin = [90;58;110;30;50;90;58;110;30;50;90;58;110;30];%火电机组下线 ru=50;rd=40;%爬坡和滑坡 %水电 Psmax = [1060;820;1243]; Psmin = [90;58;110]; rsu=0.1;rsd=0.1; cons_sgen = getConssGen(x_P_s,Psmax,Psmin,rsu,rsd,rsmax,rsmin,wwp,alfas,Horizon); cons_gen = getConsGen2(x_P_h,Phmax,Phmin,ru,rd,rgmax,rgmin,wwp,alfah,Horizon); cons = [cons, cons_gen]; cons = [cons, cons_sgen]; %风电 cons = [cons, x_P_w==muw+ww,-0.3.*muw<=ww<=0.3.*muw]; %光伏 cons = [cons, x_P_v==muv+wp,-0.3.*muv<=wp<=0.3.*muv]; % 仿射约束 x_P_hz = x_P_h-alfah.*repmat(wwp,ngen,1); x_P_sz = x_P_s-alfas.*repmat(wwp,ns,1); cons = [cons,0<= alfah <=1,0<= alfas <=1,sum(alfah)+sum(alfas) == ones(1,Horizon)];
3 程序效果
程序采用118节点系统
以上是程序优化结果图，下面为原文对照图。
以上是程序优化结果图，下面为原文对照图。4
