27 用linprog、fmincon求 解线性规划问题(matlab程序) |
您所在的位置:网站首页 › linprog函数各参数的含义 › 27 用linprog、fmincon求 解线性规划问题(matlab程序) |
1.简述
① linprog函数: 求解线性规划问题,求目标函数的最小值, [x,y]= linprog(c,A,b,Aeq,beq,lb,ub) 求最大值时,c加上负号:-c ② intlinprog函数: 求解混合整数线性规划问题, [x,y]= intlinprog(c,intcon,A,b,Aeq,beq,lb,ub) 与linprog相比,多了参数intcon,代表了整数决策变量所在的位置 优化问题中最常见的,就是线性/整数规划问题。即使赛题中有非线性目标/约束,第一想法也应是将其转化为线性。 直白点说,只要决定参加数模比赛,学会建立并求解线性/整数规划问题是非常必要的。 本期主要阐述用Matlab软件求解此类问题的一般步骤,后几期会逐步增加用Mathematica、AMPL、CPLEX、Gurobi、Python等软件求解的教程。 或许你已经听说或掌握了linprog等函数,实际上,它只是诸多求解方法中的一种,且有一定的局限性。 我的每期文章力求"阅完可上手"并"知其所以然"。因此,在讲解如何应用linprog等函数语法前,有必要先了解: 什么赛题适用线性/整数规划?如何把非线性形式线性化?如何查看某函数的语法?有哪几种求解方法?把握好这四个问题,有时候比仅仅会用linprog等函数求解更重要。 一、什么赛题适用线性/整数规划当题目中提到“怎样分配”、“XX最大/最合理”、“XX尽量多/少”等词汇时。具体有: 1. 生产安排目标:总利润最大;约束:原材料、设备限制; 2. 销售运输目标:运费等成本最低;约束:从某产地(产量有限制)运往某销地的运费不同; 3. 投资收益等目标:总收益最大;约束:不同资产配置下收益率/风险不同,总资金有限; 对于整数规划,除了通常要求变量为整数外,典型的还有指派/背包等问题(决策变量有0-1变量)。 二、如何把非线性形式线性化在比赛时,遇到非线性形式是家常便饭。此时若能够线性化该问题,绝对是你数模论文的加分项。 我在之前写的线性化文章中提到:如下非线性形式,均可实现线性化: 总的来说,具有 分段函数形式、 绝对值函数形式、 最小/大值函数形式、 逻辑或形式、 含有0-1变量的乘积形式、 混合整数形式以及 分式目标函数,均可实现 线性化。而实现线性化的主要手段主要就两点,一是引入0-1变量,二是引入很大的整数M。具体细节请参见之前写的线性化文章。 三、如何查看函数所有功能授之以鱼,不如授之以渔。 学习linprog等函数最好的方法,无疑是看Matlab官方帮助文档。本文仅是抛砖引玉地举例说明几个函数的基础用法,更多细节参见帮助文档。步骤是: 调用linprog等函数前需要事先安装“OptimizationToolbox”工具箱;在Matlab命令窗口输入“doc linprog”,便可查看语法,里面有丰富的例子;也可直接查看官方给的PDF帮助文档,后台回复“线性规划”可获取。2.代码 主程序: %% 解线性规划问题 %f(x)=-5x(1)+4x(2)+2x(3) f=[-5,4,2]; %函数系数 A=[6,-1,1;1,2,4]; %不等式系数 b=[8;10]; %不等式右边常数项 l=[-1,0,0]; %下限 u=[3,2,inf]; %上限 %%%%用linprog求解 [xol,fol]=linprog(f,A,b,[],[],l,u) %%%%用fmincon求解 x0=[0,0,0]; f1214=inline('-5*x(1)+4*x(2)+2*x(3)','x'); [xoc,foc]=fmincon(f1214,x0,A,b,[],[],l,u) 子程序: function [x,fval,exitflag,output,lambda]=linprog(f,A,B,Aeq,Beq,lb,ub,x0,options) %LINPROG Linear programming. % X = LINPROG(f,A,b) attempts to solve the linear programming problem: % % min f'*x subject to: A*x = 0. % The dual problem is % max b'*y such that A'*y + s = f, s >= 0. % % See also QUADPROG. % Copyright 1990-2018 The MathWorks, Inc. % If just 'defaults' passed in, return the default options in X % Default MaxIter, TolCon and TolFun is set to [] because its value depends % on the algorithm. defaultopt = struct( ... 'Algorithm','dual-simplex', ... 'Diagnostics','off', ... 'Display','final', ... 'LargeScale','on', ... 'MaxIter',[], ... 'MaxTime', Inf, ... 'Preprocess','basic', ... 'TolCon',[],... 'TolFun',[]); if nargin==1 && nargout 3 computeConstrViolation = true; computeFirstOrderOpt = true; % Lagrange multipliers are needed to compute first-order optimality computeLambda = true; else computeConstrViolation = false; computeFirstOrderOpt = false; computeLambda = false; end % Algorithm check: % If Algorithm is empty, it is set to its default value. algIsEmpty = ~isfield(options,'Algorithm') || isempty(options.Algorithm); if ~algIsEmpty Algorithm = optimget(options,'Algorithm',defaultopt,'fast',allDefaultOpts); OUTPUT.algorithm = Algorithm; % Make sure the algorithm choice is valid if ~any(strcmp({algIP; algDSX; algIP15b},Algorithm)) error(message('optim:linprog:InvalidAlgorithm')); end else Algorithm = algDSX; OUTPUT.algorithm = Algorithm; end % Option LargeScale = 'off' is ignored largescaleOn = strcmpi(optimget(options,'LargeScale',defaultopt,'fast',allDefaultOpts),'on'); if ~largescaleOn [linkTag, endLinkTag] = linkToAlgDefaultChangeCsh('linprog_warn_largescale'); warning(message('optim:linprog:AlgOptsConflict', Algorithm, linkTag, endLinkTag)); end % Options setup diagnostics = strcmpi(optimget(options,'Diagnostics',defaultopt,'fast',allDefaultOpts),'on'); switch optimget(options,'Display',defaultopt,'fast',allDefaultOpts) case {'final','final-detailed'} verbosity = 1; case {'off','none'} verbosity = 0; case {'iter','iter-detailed'} verbosity = 2; case {'testing'} verbosity = 3; otherwise verbosity = 1; end % Set the constraints up: defaults and check size [nineqcstr,nvarsineq] = size(A); [neqcstr,nvarseq] = size(Aeq); nvars = max([length(f),nvarsineq,nvarseq]); % In case A is empty if nvars == 0 % The problem is empty possibly due to some error in input. error(message('optim:linprog:EmptyProblem')); end if isempty(f), f=zeros(nvars,1); end if isempty(A), A=zeros(0,nvars); end if isempty(B), B=zeros(0,1); end if isempty(Aeq), Aeq=zeros(0,nvars); end if isempty(Beq), Beq=zeros(0,1); end % Set to column vectors f = f(:); B = B(:); Beq = Beq(:); if ~isequal(length(B),nineqcstr) error(message('optim:linprog:SizeMismatchRowsOfA')); elseif ~isequal(length(Beq),neqcstr) error(message('optim:linprog:SizeMismatchRowsOfAeq')); elseif ~isequal(length(f),nvarsineq) && ~isempty(A) error(message('optim:linprog:SizeMismatchColsOfA')); elseif ~isequal(length(f),nvarseq) && ~isempty(Aeq) error(message('optim:linprog:SizeMismatchColsOfAeq')); end [x0,lb,ub,msg] = checkbounds(x0,lb,ub,nvars); if ~isempty(msg) exitflag = -2; x = x0; fval = []; lambda = []; output.iterations = 0; output.constrviolation = []; output.firstorderopt = []; output.algorithm = ''; % not known at this stage output.cgiterations = []; outpussage = msg; if verbosity > 0 disp(msg) end return end if diagnostics % Do diagnostics on information so far gradflag = []; hessflag = []; constflag = false; gradconstflag = false; non_eq=0;non_ineq=0; lin_eq=size(Aeq,1); lin_ineq=size(A,1); XOUT=ones(nvars,1); funfcn{1} = []; confcn{1}=[]; diagnose('linprog',OUTPUT,gradflag,hessflag,constflag,gradconstflag,... XOUT,non_eq,non_ineq,lin_eq,lin_ineq,lb,ub,funfcn,confcn); end % Throw warning that x0 is ignored (true for all algorithms) if ~isempty(x0) && verbosity > 0 fprintf(getString(message('optim:linprog:IgnoreX0',Algorithm))); end if strcmpi(Algorithm,algIP) % Set the default values of TolFun and MaxIter for this algorithm defaultopt.TolFun = 1e-8; defaultopt.MaxIter = 85; [x,fval,lambda,exitflag,output] = lipsol(f,A,B,Aeq,Beq,lb,ub,options,defaultopt,computeLambda); elseif strcmpi(Algorithm,algDSX) || strcmpi(Algorithm,algIP15b) % Create linprog options object algoptions = optimoptions('linprog', 'Algorithm', Algorithm); % Set some algorithm specific options if isfield(options, 'InternalOptions') algoptions = setInternalOptions(algoptions, options.InternalOptions); end thisMaxIter = optimget(options,'MaxIter',defaultopt,'fast',allDefaultOpts); if strcmpi(Algorithm,algIP15b) if ischar(thisMaxIter) error(message('optim:linprog:InvalidMaxIter')); end end if strcmpi(Algorithm,algDSX) algoptions.Preprocess = optimget(options,'Preprocess',defaultopt,'fast',allDefaultOpts); algoptions.MaxTime = optimget(options,'MaxTime',defaultopt,'fast',allDefaultOpts); if ischar(thisMaxIter) && ... ~strcmpi(thisMaxIter,'10*(numberofequalities+numberofinequalities+numberofvariables)') error(message('optim:linprog:InvalidMaxIter')); end end % Set options common to dual-simplex and interior-point-r2015b algoptions.Diagnostics = optimget(options,'Diagnostics',defaultopt,'fast',allDefaultOpts); algoptions.Display = optimget(options,'Display',defaultopt,'fast',allDefaultOpts); thisTolCon = optimget(options,'TolCon',defaultopt,'fast',allDefaultOpts); if ~isempty(thisTolCon) algoptions.TolCon = thisTolCon; end thisTolFun = optimget(options,'TolFun',defaultopt,'fast',allDefaultOpts); if ~isempty(thisTolFun) algoptions.TolFun = thisTolFun; end if ~isempty(thisMaxIter) && ~ischar(thisMaxIter) % At this point, thisMaxIter is either % * a double that we can set in the options object or % * the default string, which we do not have to set as algoptions % is constructed with MaxIter at its default value algoptions.MaxIter = thisMaxIter; end % Create a problem structure. Individually creating each field is quicker % than one call to struct problem.f = f; problem.Aineq = A; problem.bineq = B; problem.Aeq = Aeq; problem.beq = Beq; problem.lb = lb; problem.ub = ub; problem.options = algoptions; problem.solver = 'linprog'; % Create the algorithm from the options. algorithm = createAlgorithm(problem.options); % Check that we can run the problem. try problem = checkRun(algorithm, problem, 'linprog'); catch ME throw(ME); end % Run the algorithm [x, fval, exitflag, output, lambda] = run(algorithm, problem); % If exitflag is {NaN, }, this means an internal error has been % thrown. The internal exit code is held in exitflag{2}. if iscell(exitflag) && isnan(exitflag{1}) handleInternalError(exitflag{2}, 'linprog'); end end output.algorithm = Algorithm; % Compute constraint violation when x is not empty (interior-point/simplex presolve % can return empty x). if computeConstrViolation && ~isempty(x) output.constrviolation = max([0; norm(Aeq*x-Beq, inf); (lb-x); (x-ub); (A*x-B)]); else output.constrviolation = []; end % Compute first order optimality if needed. This information does not come % from either qpsub, lipsol, or simplex. if exitflag ~= -9 && computeFirstOrderOpt && ~isempty(lambda) output.firstorderopt = computeKKTErrorForQPLP([],f,A,B,Aeq,Beq,lb,ub,lambda,x); else output.firstorderopt = []; end 3.运行结果
|
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |