电力系统中新型预测双二元变量机组组合问题（Matlab代码实现）

目录

1 概述

2 参考文献

3 Matlab代码实现

---

1 概述

高效求解大规模 SCUC 问题的关键在于削减其规模。文献[1]表明，安全约束机组组合问题中

大量的故障态安全约束是冗余且无效的，不会对SCUC 问题的最优解产生影响。因此，可以通过辨

识、删除这些冗余约束，将 SCUC 问题的规模进行大幅度削减，从而缩短 SCUC 问题的求解时间，以实现 SCUC 问题的快速求解。文献[19]基于一系列小规模 MILP 问题的解析解给出了辨识冗余安全约束的充分不必要条件，但是它所研究的 SCUC 模型忽略了故障态安全约束，不够全面。文献[2]提出通过简化故障态集合来削减故障态安全约束，但是经过简化后的故障态集合仍然包含很多冗余的故障态安全约束，其作用较为有限。文献提出了一 种基于迭代过程约束削减方法，它逐渐将未满足的约束条件添加到 SCUC 模型中，并进行迭代求解。 但是，该方法在某些情况下需要进行多次迭代，反而会延长求解时间。部分代码：

function [ miqp_i, miqp ] = miqpUC_no_proj( dataUC )

N = dataUC.N;   T = dataUC.T;
x_L = cell(N,1);    x_U = cell(N,1);
A = cell(N,1);      b = cell(N,1);
B = cell(N,1);      B_wan = cell(N,1);
Q_UC = cell(N,1);   c_UC = cell(N,1);
ctype = cell(N,1);

DiagE_T_T = sparse(1:T,1:T,1);
for i = 1:N   % 按照机组分块形成约束，对机组i，生成所有其相关约束作为一块
    x_L{i} = sparse(2*T,1);
    x_U{i} = [ones(1*T,1); inf* ones(T,1)];
    
    % unit generation limits 
    A_low_up = [ -sparse(1:T,1:T,dataUC.p_up(i)), DiagE_T_T;
                 sparse(1:T,1:T,dataUC.p_low(i)), -DiagE_T_T;
                ];
    b_low_up =  sparse(2*T,1);   
    
    
    
    % form A_i  b_i
    A{i} = [    A_low_up;
        ];
    b{i} = [    b_low_up;
        ];
    
    % 对功率平衡和旋转备用约束形成约束
    B{i} =  [ sparse(1:T,1:T,-dataUC.p_up(i)), sparse(T, 1*T)];
    B_wan{i} = [sparse(T,T), DiagE_T_T];

    % 目标函数
    c_UC{i} = [ dataUC.alpha(i) * ones(T,1); 
                dataUC.beta(i) * ones(T,1); 
                ];
    Q_UC{i} = blkdiag(sparse(T,T), sparse(1:T,1:T, 2*dataUC.gamma(i)));
    
    ctype_ascii_miqp = [66*ones(1,T), 67*ones(1,T)];
    ctype{i} = char(ctype_ascii_miqp);    
end

% 形成功率平衡和备用约束的右端项
c = -(dataUC.PD + dataUC.spin);
c_wan = dataUC.PD;

% 返回分块的模型数据
miqp_i.A = A;           miqp_i.b = b;
miqp_i.B = B;           miqp_i.c = c;
miqp_i.B_wan = B_wan;   miqp_i.c_wan = c_wan;
miqp_i.c_UC = c_UC;     miqp_i.Q_UC = Q_UC;
miqp_i.x_L = x_L;       miqp_i.x_U = x_U;
miqp_i.N = N;           miqp_i.T = T;
miqp_i.ctype = ctype;

% 形成并返回 完整MIQP模型的参数
miqp.A = [];           miqp.b = [];
miqp.B = [];           miqp.c = c;
miqp.B_wan = [];   miqp.c_wan = c_wan;
miqp.c_UC = [];     miqp.Q_UC = [];
miqp.x_L = [];       miqp.x_U = [];
miqp.ctype = [];
for i = 1:N
   miqp.A = blkdiag(miqp.A, A{i});
   miqp.b = [miqp.b; b{i}];
   miqp.B = [miqp.B, B{i}];
   miqp.B_wan = [miqp.B_wan, B_wan{i}];
   miqp.c_UC = [miqp.c_UC; c_UC{i}];
   miqp.Q_UC = blkdiag(miqp.Q_UC, Q_UC{i});
   miqp.x_L = [miqp.x_L; x_L{i}];
   miqp.x_U = [miqp.x_U; x_U{i}];
   miqp.ctype = [miqp.ctype, ctype{i}];
end
miqp.A_all =  [     miqp.A;
                    miqp.B;
%                     miqp.B_wan;
    ];
miqp.lhs_all = [    -inf * ones(size(miqp.A,1),1);
                    -inf * ones(size(miqp.B,1),1);
%                     miqp.c_wan;
    ];
miqp.rhs_all = [    miqp.b;
                    miqp.c;
%                     miqp.c_wan;
    ];
miqp.N = N;
miqp.T = T;

%%%% end of function
end

function [ miqp_i, miqp ] = miqpUC_no_proj( dataUC )

N = dataUC.N;   T = dataUC.T;
x_L = cell(N,1);    x_U = cell(N,1);
A = cell(N,1);      b = cell(N,1);
B = cell(N,1);      B_wan = cell(N,1);
Q_UC = cell(N,1);   c_UC = cell(N,1);
ctype = cell(N,1);

DiagE_T_T = sparse(1:T,1:T,1);
for i = 1:N   % 按照机组分块形成约束，对机组i，生成所有其相关约束作为一块
    x_L{i} = sparse(2*T,1);
    x_U{i} = [ones(1*T,1); inf* ones(T,1)];
    
    % unit generation limits 
    A_low_up = [ -sparse(1:T,1:T,dataUC.p_up(i)), DiagE_T_T;
                 sparse(1:T,1:T,dataUC.p_low(i)), -DiagE_T_T;
                ];
    b_low_up =  sparse(2*T,1);   
    
    
    
    % form A_i  b_i
    A{i} = [    A_low_up;
        ];
    b{i} = [    b_low_up;
        ];
    
    % 对功率平衡和旋转备用约束形成约束
    B{i} =  [ sparse(1:T,1:T,-dataUC.p_up(i)), sparse(T, 1*T)];
    B_wan{i} = [sparse(T,T), DiagE_T_T];

    % 目标函数
    c_UC{i} = [ dataUC.alpha(i) * ones(T,1); 
                dataUC.beta(i) * ones(T,1); 
                ];
    Q_UC{i} = blkdiag(sparse(T,T), sparse(1:T,1:T, 2*dataUC.gamma(i)));
    
    ctype_ascii_miqp = [66*ones(1,T), 67*ones(1,T)];
    ctype{i} = char(ctype_ascii_miqp);    
end

% 形成功率平衡和备用约束的右端项
c = -(dataUC.PD + dataUC.spin);
c_wan = dataUC.PD;

% 返回分块的模型数据
miqp_i.A = A;           miqp_i.b = b;
miqp_i.B = B;           miqp_i.c = c;
miqp_i.B_wan = B_wan;   miqp_i.c_wan = c_wan;
miqp_i.c_UC = c_UC;     miqp_i.Q_UC = Q_UC;
miqp_i.x_L = x_L;       miqp_i.x_U = x_U;
miqp_i.N = N;           miqp_i.T = T;
miqp_i.ctype = ctype;

% 形成并返回 完整MIQP模型的参数
miqp.A = [];           miqp.b = [];
miqp.B = [];           miqp.c = c;
miqp.B_wan = [];   miqp.c_wan = c_wan;
miqp.c_UC = [];     miqp.Q_UC = [];
miqp.x_L = [];       miqp.x_U = [];
miqp.ctype = [];
for i = 1:N
   miqp.A = blkdiag(miqp.A, A{i});
   miqp.b = [miqp.b; b{i}];
   miqp.B = [miqp.B, B{i}];
   miqp.B_wan = [miqp.B_wan, B_wan{i}];
   miqp.c_UC = [miqp.c_UC; c_UC{i}];
   miqp.Q_UC = blkdiag(miqp.Q_UC, Q_UC{i});
   miqp.x_L = [miqp.x_L; x_L{i}];
   miqp.x_U = [miqp.x_U; x_U{i}];
   miqp.ctype = [miqp.ctype, ctype{i}];
end
miqp.A_all =  [     miqp.A;
                    miqp.B;
%                     miqp.B_wan;
    ];
miqp.lhs_all = [    -inf * ones(size(miqp.A,1),1);
                    -inf * ones(size(miqp.B,1),1);
%                     miqp.c_wan;
    ];
miqp.rhs_all = [    miqp.b;
                    miqp.c;
%                     miqp.c_wan;
    ];
miqp.N = N;
miqp.T = T;

%%%% end of function
end

2 参考文献

[1]李建钊,谢敏,李舒佳,林盛振,黄彬彬.考虑CVaR的机组组合和多场景备用决策联合优化[J].南方能源建设,2021,8(04):50-65.DOI:10.16516/j.gedi.issn2095-8676.2021.04.008.

[2]石志伟.计及FACTS校正的N-k故障约束鲁棒机组组合[J].南方电网技术,2021,15(12):11-19.DOI:10.13648/j.cnki.issn1674-0629.2021.12.002.

3 Matlab代码实现