27 用linprog、fmincon求 解线性规划问题(matlab程序)

您所在的位置:网站首页 linprog函数各参数的含义 27 用linprog、fmincon求 解线性规划问题(matlab程序)

27 用linprog、fmincon求 解线性规划问题(matlab程序)

2024-01-15 17:17| 来源: 网络整理| 查看: 265

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