(文章复现)微电网两阶段鲁棒优化经济调度方法(含matlab代码)

一、文献中存在的笔误

文献中两阶段鲁棒优化的模型的构建是基于下面这个确定性优化：

但是我复现的时候发现，这个模型其实是有点小问题的，第三行约束条件写的是Ky=0，但是对应的约束条件中等式右边并不全为零，所以按照常理，这一个约束应该写成Ky=k。但不确定这个文章是做了线性变换还是笔误。我也在这里纠结了好久。

后面不知怎么的，我想着给文章通信作者发邮件询问一下，没想到真的收到回复了。刘一欣老师说这里确实是个笔误，等式右边应还有一个列向量。特别感谢刘老师可以回复我这个小菜鸡，不然我要在这上面纠结很久。

我询问的邮件：

刘一欣老师的回复：

我邮件中提到的另一个问题，是关于这个对偶问题的：

因为根据强对偶定理，对偶变量对应的约束为等式约束，所以化为对偶形式时，对偶变量不应该含有符号约束，也就是≥0这个约束要删去。

同时因为约束Ky=0这一项要改成Ky=k，所以子问题的目标函数也要加上这一项。

解决了上面这些问题，后面的建模就相对没有那么困难了。

二、决策变量和约束矩阵

决策变量x是一个48维的列向量：

决策变量y是一个240维的列向量： 

 c和y一样都是一个240维的列向量：

 D是一个144行240列的矩阵，d是一个144维的列向量，Dy≥d可以写成矩阵形式：

表示24阶的单位矩阵，表示下三角和对角线全为1的矩阵，上三角全为0的矩阵。

K是一个50行240列的矩阵，k是一个50维的列向量，Ky=k可以写成矩阵形式：

式中，表示元素全为1的24维行向量。

F是一个96行48列的矩阵，G是一个96行240列的矩阵，h为96维的列向量，Fx+Gy≥h可以写成矩阵形式：

 

这样就实现了完成了所有系数矩阵的输入。

其他没有在矩阵中涉及的约束条件，用y>0表示即可。

三、matlab代码

%% 清理内存空间
clc
clear
close all
%% 随机生成初始的光伏出力和负荷大小
u_sharp=creat_situation;
%% 迭代求解两阶段鲁棒优化问题
it_max=10;%最大迭代次数
LB=zeros(1,it_max);
UB=zeros(1,it_max);
e=zeros(1,it_max);
for it=1:it_max
    [x,LB(it),~]=master_problem(u_sharp);
    [u_sharp,UB(it)]=sub_problem(x);
    e(it)=UB(it)-LB(it);
    if abs(e(it))<0.01
        break
    end
end
%% 画图
index=[1 1 1 1 1 1 1];
%共有七个图,1-算法迭代情况,2-光伏出力,3-负荷大小,4-微型燃气轮机出力,5-微网购售电,6-储能充放电,7-需求响应负荷
%index对应位置为1表示画图，为0表示跳过
plot_figure(index,u_sharp,e,it_max);

上面只给出了主函数，完整代码戳这里：

原创代码：基于matlab+yalmip+cplex实现的微电网两阶段鲁棒优化

       

四、运行结果

 运行结果一共有七个图，效果如下：

 

 

 

因为没法获取原文献的数据，所以具体结果和文献不会完全一致，但是原理都是一样的。 

PS：统一回复，博客和代码发布后，有很多朋友在私信和我交流，原文使用强对偶定理求解的模型存在一些问题，参数不同，求解器版本不同等情况可能使代码没办法收敛。我调试之后，改用KKT条件改写了子问题，可以保证不同参数、不同求解器版本下模型的收敛性。