混合整数规划的机组组合附Matlab代码

✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知。

更多Matlab仿真内容点击👇

智能优化算法神经网络预测雷达通信无线传感器

信号处理图像处理路径规划元胞自动机无人机电力系统

⛄ 内容介绍

机组组合问题要求基于已知的系统数据，求解计划时间内机组决策变量的最优组合，使得系统总成本达到最小。该问题的决策变量由两类，一类是各时段机组的启停状态，为整数变量，0表示关停，1表示启动；另一类是各时段机组的出力，为连续变量。

机组组合问题属于规划问题，即要在决策变量的可行解空间里找到一组最优解，使得目标函数尽可能取得极值。对于混合整数规划，常用的方法有分支定界法，benders分解等。CPLEX提供了快速的MIP求解方法，对于数学模型已知的问题，只需要按照程序规范在MATLAB中编写程序化模型，调用CPLEX求解器，即可进行求解。

下文介绍机组组合优化的数学模型。

校验程序的算例基于IEEE-30节点标准测试系统，系统接线图如图1。系统包含30个节点，6台发电机组。要求确定系统最优机组组合，使得系统各机组总运行成本（煤耗成本+启停成本）最小化。

图1. IEEE-30节点测试系统接线

已知：给定系统数据包括如下：（见附件testsystem.xls）

1)线路网络参数

2)机组参数

3)各节点各时段负荷曲线（24小时）

注意：附件中的数据均基于标幺化系统得到，因此电力电量参数、网络参数等都为标幺值，无量纲。还要注意附件中煤耗系数a,b,c的单位为吨，因此计算煤耗成本还需换算为价格，设燃煤价格为100$/吨。

求解：机组组合结果，即机组各时段启停计划、机组各时段最优出力，以及内含的各时段的直流潮流等。

⛄ 部分代码

function diagnostic = solvesdp(varargin)

%SOLVESDP Computes solution to optimization problem

% DIAGNOSTIC = SOLVESDP(F,h,options) is the common command to

% solve optimization problems of the following kind

% min h

% subject to

% F >=(<=,==) 0

% NOTES

% Despite the name, SOLVESDP is the interface for solving all

% supported problem classes (LP, QP, SOCP, SDP, BMI, MILP, MIQP,...)

% To obtain solution for a variable, use DOUBLE.

% To obtain dual variable for a constraint, use DUAL.

% See YALMIPERROR for error codes returned in output.

% OUTPUT

% diagnostic : Diagnostic information

% INPUT

% F : Object describing the constraints. Can be [].

% h : SDPVAR object describing the objective h(x). Can be [].

% options : Options structure. See SDPSETTINGS. Can be [].

% EXAMPLE

% A = randn(15,5);b = rand(15,1)*5;c = randn(5,1);

% x = sdpvar(5,1);

% solvesdp([x>=0, A*x<=b],c'*x);double(x)

% See also DUAL, @SDPVAR/DOUBLE, SDPSETTINGS, YALMIPERROR

yalmiptime = clock; % Let us see how much time we spend

% Avoid warning

if length(varargin)>=2

if isa(varargin{2},'double')

varargin{2} = [];

end

if length(varargin)>=2

if isa(varargin{2},'sdpvar') && prod(size(varargin{2}))>1

% Several objectives

diagnostic = solvesdp_multiple(varargin{:});

return

end

% *********************************

% CHECK INPUT

% *********************************

nargin = length(varargin);

if nargin<1

help solvesdp

return

else

F = varargin{1};

if isa(F,'constraint')

F = lmi(F);

end

if isa(F,'lmi')

F = flatten(F);

end

if isa(F,'sdpvar')

% We do allow sloppy coding of logic constraints, i.e writing a

% constraints as [a|b true(a)]

Fnew = [];

for i = 1:length(F)

if length(getvariables(F(i)))>1

Fnew = nan;

break

end

operator = yalmip('extstruct',getvariables(F(i)));

if isempty(operator)

Fnew = nan;

break

end

if length(operator)>1

Fnew = nan;

break

end

if ~strcmp(operator.fcn,'or')

Fnew = nan;

break

end

Fnew = Fnew + (true(F(i)));

end

if isnan(Fnew)

error('First argument (F) should be a constraint object.');

else

F = Fnew;

end

elseif isempty(F)

F = lmi([]);

elseif ~isa(F,'lmi')

error('First argument (F) should be a constraint object.');

end

if nargin>=2

h = varargin{2};

if isa(h,'double')

h = [];

end

if ~(isempty(h) | isa(h,'sdpvar') | isa(h,'logdet') | isa(h,'ncvar'))

error('Second argument (the objective function h) should be an sdpvar or logdet object (or empty).');

end

if isa(h,'logdet')

logdetStruct.P = getP(h);

logdetStruct.gain = getgain(h);

h = getcx(h);

if isempty(F)

F = ([]);

end

else

logdetStruct = [];

end

else

logdetStruct = [];

h = [];

end

if ~isempty(F)

if any(is(F,'sos'))

diagnostic = solvesos(varargin{:});

return

end

if isa(h,'sdpvar')

if is(h,'complex')

error('Complex valued objective does not make sense.');

end

if nargin>=3

options = varargin{3};

if ~(isempty(options) | isa(options,'struct'))

error('Third argument (options) should be an sdpsettings struct (or empty).');

end

if isempty(options)

options = sdpsettings;

end

else

options = sdpsettings;

end

options.solver = lower(options.solver);

% If user has logdet term, but no preference on solver, we try to hook up

% with SDPT3 if possible.

if ~isempty(logdetStruct)

if strcmp(options.solver,'')

% options.solver = 'sdpt3,*';

end

% Call chance solver?

if length(F) > 0

rand_declarations = is(F,'random');

if any(rand_declarations)

% diagnostic = solverandom(F(find(~rand_declarations)),h,options,recover(getvariables(sdpvar(F(find(unc_declarations))))));

return

end

% Call robust solver?

if length(F) > 0

unc_declarations = is(F,'uncertain');

if any(unc_declarations)

diagnostic = solverobust(F(find(~unc_declarations)),h,options,recover(getvariables(sdpvar(F(find(unc_declarations))))));

return

end

if isequal(options.solver,'mpt') | nargin>=4

solving_parametric = 1;

else

solving_parametric = 0;

end

% Just for safety

if isempty(F) & isempty(logdetStruct)

F = lmi;

end

if any(is(F,'sos'))

error('You have SOS constraints. Perhaps you meant to call SOLVESOS.');

end

% Super stupido

if length(F) == 0 & isempty(h) & isempty(logdetStruct)

diagnostic.yalmiptime = 0;

diagnostic.solvertime = 0;

diagnostic.info = 'No problems detected (YALMIP)';

diagnostic.problem = 0;

diagnostic.dimacs = [NaN NaN NaN NaN NaN NaN];

return

end

% Dualize the problem?

if ~isempty(F)

if options.dualize == -1

sdp = find(is(F,'sdp'));

if ~isempty(sdp)

if all(is(F(sdp),'sdpcone'))

options.dualize = 1;

end

if options.dualize == 1

[Fd,objd,aux1,aux2,aux3,complexInfo] = dualize(F,h,[],[],[],options);

options.dualize = 0;

diagnostic = solvesdp(Fd,-objd,options);

if ~isempty(complexInfo)

for i = 1:length(complexInfo.replaced)

n = size(complexInfo.replaced{i},1);

re = 2*double(complexInfo.new{i}(1:n,1:n));

im = 2*double(complexInfo.new{i}(1:n,n+1:end));

im=triu((im-im')/2)-(triu((im-im')/2))';

assign(complexInfo.replaced{i},re + sqrt(-1)*im);

end

return

end

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

% DID WE SELECT THE MOMENT SOLVER

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

if isequal(options.solver,'moment')

if ~isempty(logdetStruct)

error('Cannot dualize problems with logaritmic objective')

end

options.solver = options.moment.solver;

[diagnostic,x,momentdata] = solvemoment(F,h,options,options.moment.order);

diagnostic.momentdata = momentdata;

diagnostic.xoptimal = x;

return

end

% ******************************************

% COMPILE IN GENERALIZED YALMIP FORMAT

% ******************************************

[interfacedata,recoverdata,solver,diagnostic,F,Fremoved,ForiginalQuadratics] = compileinterfacedata(F,[],logdetStruct,h,options,0,solving_parametric);

% ******************************************

% FAILURE?

% ******************************************

if ~isempty(diagnostic)

diagnostic.yalmiptime = etime(clock,yalmiptime);

return

end

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

% DID WE SELECT THE LMILAB SOLVER WITH A KYP

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

if strcmpi(solver.tag,'lmilab') & any(is(F,'kyp'))

[diagnostic,failed] = calllmilabstructure(F,h,options);

if ~failed % Did this problem pass (otherwise solve using unstructured call)

diagnostic.yalmiptime = etime(clock,yalmiptime)-diagnostic.solvertime;

return

end

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

% DID WE SELECT THE KYPD SOLVER

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

if strcmpi(solver.tag,'kypd')

diagnostic = callkypd(F,h,options);

diagnostic.yalmiptime = etime(clock,yalmiptime)-diagnostic.solvertime;

return

end

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

% DID WE SELECT THE STRUL SOLVER

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

if strfind(solver.tag,'STRUL')

diagnostic = callstrul(F,h,options);

diagnostic.yalmiptime = etime(clock,yalmiptime)-diagnostic.solvertime;

return

end

% ******************************************

% DID WE SELECT THE BMILIN SOLVER (obsolete)

% ******************************************

if strcmpi(solver.tag,'bmilin')

diagnostic = callbmilin(F,h,options);

return

end

% ******************************************

% DID WE SELECT THE BMIALT SOLVER (obsolete)

% ******************************************

if strcmp(solver.tag,'bmialt')

diagnostic = callbmialt(F,h,options);

return

end

%******************************************

% DID WE SELECT THE MPT solver (backwards comb)

%******************************************

actually_save_output = interfacedata.options.savesolveroutput;

if strcmpi(solver.tag,'mpt-2') | strcmpi(solver.tag,'mpt-3') | strcmpi(solver.tag,'mpcvx') | strcmpi(solver.tag,'mplcp')

interfacedata.options.savesolveroutput = 1;

if isempty(interfacedata.parametric_variables)

if (nargin < 4 | ~isa(varargin{4},'sdpvar'))

error('You must specify parametric variables.')

else

interfacedata.parametric_variables = [];

for i = 1:length(varargin{4})

interfacedata.parametric_variables = [interfacedata.parametric_variables;find(ismember(recoverdata.used_variables,getvariables(varargin{4}(i))))];

end

if isempty(varargin{5})

interfacedata.requested_variables = [];

else

interfacedata.requested_variables = [];

for i = 1:length(varargin{5})

interfacedata.requested_variables = [interfacedata.requested_variables;find(ismember(recoverdata.used_variables,getvariables(varargin{5}(i))))];

end

% *************************************************************************

% Just return the YALMIP model. Used when solving multiple objectives

% *************************************************************************

if isfield(options,'pureexport')

interfacedata.recoverdata = recoverdata;

diagnostic = interfacedata;

return

end

% *************************************************************************

% TRY TO SOLVE PROBLEM

% *************************************************************************

if options.debug

eval(['output = ' solver.call '(interfacedata);']);

else

try

eval(['output = ' solver.call '(interfacedata);']);

catch

output.Primal = zeros(length(interfacedata.c),1)+NaN;

output.Dual = [];

output.Slack = [];

output.solvertime = nan;

output.solverinput = [];

output.solveroutput = [];

output.problem = 9;

output.infostr = yalmiperror(output.problem,lasterr);

end

if options.dimacs

try

b = -interfacedata.c;

c = interfacedata.F_struc(:,1);

A = -interfacedata.F_struc(:,2:end)';

x = output.Dual;

y = output.Primal;

% FIX this nonlinear crap (return variable type in

% compileinterfacedata)

if options.relax == 0 & any(full(sum(interfacedata.monomtable,2)~=0))

if ~isempty(find(sum(interfacedata.monomtable | interfacedata.monomtable,2)>1))

z=real(exp(interfacedata.monomtable*log(y+eps)));

y = z;

end

if isfield(output,'Slack')

s = output.Slack;

else

s = [];

end

dimacs = computedimacs(b,c,A,x,y,s,interfacedata.K);

catch

dimacs = [nan nan nan nan nan nan];

end

else

dimacs = [nan nan nan nan nan nan];

end

% ********************************

% ORIGINAL COORDINATES

% ********************************

output.Primal = recoverdata.x_equ+recoverdata.H*output.Primal;

% ********************************

% OUTPUT

% ********************************

diagnostic.yalmiptime = etime(clock,yalmiptime)-output.solvertime;

diagnostic.solvertime = output.solvertime;

try

diagnostic.info = output.infostr;

catch

diagnostic.info = yalmiperror(output.problem,solver.tag);

end

diagnostic.problem = output.problem;

if options.dimacs

diagnostic.dimacs = dimacs;

end

% Some more info is saved internally

solution_internal = diagnostic;

solution_internal.variables = recoverdata.used_variables(:);

solution_internal.optvar = output.Primal;

if ~isempty(interfacedata.parametric_variables)

diagnostic.mpsol = output.solveroutput;

options.savesolveroutput = actually_save_output;

end;

if interfacedata.options.savesolveroutput

diagnostic.solveroutput = output.solveroutput;

end

if interfacedata.options.savesolverinput

diagnostic.solverinput = output.solverinput;

end

if interfacedata.options.saveyalmipmodel

diagnostic.yalmipmodel = interfacedata;

end

if options.warning & warningon & isempty(findstr(diagnostic.info,'No problems detected'))

disp(['Warning: ' output.infostr]);

end

if ismember(output.problem,options.beeponproblem)

try

beep; % does not exist on all ML versions

catch

end

% And we are done! Save the result

if ~isempty(output.Primal)

if size(output.Primal,2)>1

for j = 1:size(output.Primal,2)

temp = solution_internal;

temp.optvar = temp.optvar(:,j);

yalmip('setsolution',temp,j);

end

else

yalmip('setsolution',solution_internal);

end

if interfacedata.options.saveduals & solver.dual

if isempty(interfacedata.Fremoved) | (nnz(interfacedata.Q)>0)

try

setduals(F,output.Dual,interfacedata.K);

catch

end

else

try

% Duals related to equality constraints/free variables

% have to be recovered b-A*x-Ht == 0

b = -interfacedata.oldc;

A = -interfacedata.oldF_struc(1+interfacedata.oldK.f:end,2:end)';

H = -interfacedata.oldF_struc(1:interfacedata.oldK.f,2:end)';

x = output.Dual;

b_equ = b-A*x;

newdual = H\b_equ;

setduals(interfacedata.Fremoved + F,[newdual;output.Dual],interfacedata.oldK);

catch

% this is a new feature...

disp('Dual recovery failed. Please report this issue.');

end

% Hack to recover original QCQP duals from gurobi

if strcmp(solver.tag,'GUROBI-GUROBI')

if length(ForiginalQuadratics) > 0

if isfield(output,'qcDual')

if length(output.qcDual) == length(ForiginalQuadratics)

Ktemp.l = length(output.qcDual);

Ktemp.f = 0;

Ktemp.q = 0;

Ktemp.s = 0;

Ktemp.r = 0;

setduals(ForiginalQuadratics,-output.qcDual,Ktemp);

end

function yesno = warningon

s = warning;

yesno = isequal(s,'on');

⛄ 运行结果

⛄ 参考文献

[1]程杉王贤宁冯毅煁王睿娟. 基于CPLEX与MATLAB的电动汽车充电站优化调度仿真系统[J]. 电网与清洁能源, 2018, 034(001):123-127,136.

❤️ 关注我领取海量matlab电子书和数学建模资料

❤️部分理论引用网络文献，若有侵权联系博主删除