MATLAB课件第十一章线性极值
第十一章 线性极值MATLAB提供了很多求极值（或最优值）的命令函数，既可以求无条件的极值，也可求有条件的极值，其中，条件可以是不等式，也可以是等式的，可以是线性的，也可以是非线性的，甚至可以是多个条件，目标函数可以是线性的，也可以是非线性的，总之，MATLAB针对不同的类型，采用不同的函数命令去求解，以下将分类型来做些简单的介绍1线性极值（又称线性规划）1.1线性规划模型规划问题研究的对象大体可以分为两大类：一类是在现有的人、财、物等资源的条件下，研究如何合理的计划、安排，可使得某一目标达到最大，如产量、利润目标等；另一类是在任务确定后，如何计划、安排，使用最低限度的人、财等资源，去实现该任务，如使成本、费用最小等这两类问题从本质上说是相同的，即都在一组约束条件下，去实现某一个目标的最优（最大或最小）线性规划研究的问题要求目标与约束条件函数都是线性的，而目标函数只能是一个在经济管理问题中，大量问题是线性的，有的也可以转化为线性的，从而使线性规划有极大的应用价值线性规划模型包含3个要素：（1）决策变量. 问题中需要求解的那些未知量，一般用n维向量表示2）目标函数. 通常是问题需要优化的那个目标的数学表达式，它是决策变量x的线性函数。
3）约束条件. 对决策变量的限制条件，即x的允许取值范围，它通常是x的一组线性不等式或线性等式线性规划问题的数学模型一般可表示为： min（max） f T X s.t A X≤b Aeq X =beqlb≤X≤ub 其中X为n维未知向量，f T=[f1,f2,…fn]为目标函数系数向量，小于等于约束系数矩阵A为m×n矩阵，b为其右端m维列向量，Aeq为等式约束系数矩阵，beq为等式约束右端常数列向量lb,ub为自变量取值上界与下界约束的n维常数向量特别注意：当我们用MATLAB软件作优化问题时，所有求maxf 的问题化为求min(-f )来作约束g i (x)≥0，化为 –g i≤0来做1.2．线性规划问题求最优解函数： 调用格式： x=linprog(f,A,b) x=linprog(f,A,b,Aeq,beq) x=linprog(f,A,b,Aeq,beq,lb,ub) x=linprog(f,A,b,Aeq,beq,lb,ub,x0) x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options) [x,fval]=linprog(…) [x, fval, exitflag]=linprog(…) [x, fval, exitflag, output]=linprog(…) [x, fval, exitflag, output, lambda]=linprog(…) 说明：x=linprog(f,A,b)返回值x为最优解向量。
x=linprog(f,A,b,Aeq,beq) 作有等式约束的问题若没有不等式约束，则令A=[ ]、b=[ ] x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options) 中lb ,ub为变量x的下界和上界，x0为初值点，options为指定优化参数进行最小化 [x,fval]=linprog(…) 左端 fval 返回解x处的目标函数值[x,fval,exitflag,output,lambda]=linprog(f,A,b, Aeq,beq,lb,ub,x0) 的输出部分： exitflag 描述函数计算的退出条件：若为正值，表示目标函数收敛于解x处；若为负值，表示目标函数不收敛；若为零值，表示已经达到函数评价或迭代的最大次数Output为关于优化的一些信息Lambda为解x的Lagrange乘子例1】求解线性规划问题： max f=2x1+5x2 s.t 先将目标函数转化成最小值问题：min(-f)=- 2x1-5x2具体程序如下：f=[-2 -5];A=[1 0;0 1;1 2];b=[4;3;8];lb=[0 0];[x,fval]=linprog(f,A,b,[],[],lb)f=fval*(-1)运行结果： x = 2 3 fval = -19.0000maxf = 19【例2】：minf=5x1-x2+2x3+3x4-8x5s.t –2x1+x2-x3+x4-3x5≤6 2x1+x2-x3+4x4+x5≤7 0≤xj≤15 j=1,2,3,4,5编写以下程序：f=[5 -1 2 3 -8];A=[-2 1 -1 1 -3;2 1 -1 4 1];b=[6;7];lb=[0 0 0 0 0];ub=[15 15 15 15 15];[x,fval]=linprog(f,A,b,[],[],lb,ub) 运行结果：x = 0.0000 0.0000 8.0000 0.0000 15.0000minf = -104【例3】：假设某厂计划生产甲、乙两种产品，现库存主要材料有A类3600公斤，B类2000公斤，C类3000公斤。
每件甲产品需用材料A类9公斤，B类4公斤，C类3公斤每件乙产品，需用材料A类4公斤，B类5公斤，C类10公斤甲单位产品的利润70元，乙单位产品的利润120元问如何安排生产，才能使该厂所获的利润最大建立数学模型：设x1、x2分别为生产甲、乙产品的件数f为该厂所获总润 max f=70x1+120x2 s.t 9x1+4x2≤3600 4x1+5x2≤2000 3x1+10x2≤3000 x1,x2≥0将其转换为标准形式：min f=-70x1-120x2 s.t 9x1+4x2≤3600 4x1+5x2≤2000 3x1+10x2≤3000 x1,x2≥0 编写以下程序： f=[-70 -120]; A=[9 4 ;4 5;3 10 ]; b=[3600;2000;3000]; lb=[0 0]; [x,fval,exitflag]=linprog(f,A,b,[],[],lb); x, exitflag,maxf=-fval 运行结果： x = 200.0000 240.0000exitflag = 1maxf = 4.2800e+004【例4】：某公司有一批资金用于4个工程项目的投资，其投资各项目时所得的净收益(投入资金百分比)如下表：工程项目收益表工程项目ABCD收益(%)1510812由于某种原因，决定用于项目A的投资不大于其他各项投资之和而用于项目B和C的投资要大于项目D的投资。
试确定该公司收益最大的投资分配方案建立数学模型：设x1、 x2 、x3 、x4分别代表用于项目A、B、C、D的投资百分数 max f=0.15x1+0.1x2+0.08 x3+0.12 x4 s.t x1-x2- x3- x4≤0 x2+ x3- x4≥0 x1+x2+x3+ x4=1 xj≥0 j=1,2,3,4将其转换为标准形式：min z=-0.15x1-0.1x2-0.08 x3-0.12 x4 s.t x1-x2- x3- x4≤0 -x2- x3+ x4≤0 x1+x2+x3+ x4=1 xj≥0 j=1,2,3,4 编写程序： f = [-0.15;-0.1;-0.08;-0.12];A = [1 -1 -1 -1;0 -1 -1 1];b = [0; 0];Aeq=[1 1 1 1];beq=[1];lb = zeros(4,1); [x,fval,exitflag] = linprog(f,A,b,Aeq,beq,lb) fmax=-fval 运行结果：x = 0.5000 0.2500 0.0000 0.2500fval = -0.1300exitflag = 1fmax =0.1300 即4个项目的投资百分数分别为50%，25%，0, 25%时可使该公司获得最大的收益，其最大收益可到达13%。
过程正常收敛例5】：有A、B、C三个食品加工厂，负责供给甲、乙、丙、丁四个市场三个厂每天生产食品箱数上限如下表：工厂ABC生产数604050四个市场每天的需求量如下表：市场甲乙丙丁需求量20353334从各厂运到各市场的运输费(元/每箱)由下表给出：收点发点市 场甲乙丙丁工厂A2132B1321C3411求在基本满足供需平衡的约束条件下使总运输费用最小建立数学模型：设ai j为由工厂i运到市场j的费用，xi j 是由工厂i运到市场j的箱数bi是工厂i的产量，dj是市场j的需求量 b= ( 60 40 50 ) d= ( 20 35 33 34 ) s.t x i j≥0 编写程序： AA=[2 1 3 2;1 3 2 1;3 4 1 1]; f=AA(:); A=[ 1 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 1]; Aeq=[1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1]; b=[60;40;50]; beq=[20;35;33;34]; lb=zeros(12,1); [x,fval,exitflag]=linprog(f,A,b,Aeq,beq,lb) 运行结果： x = 0.0000 20.0000 0.0000 35.0000 0.0000 0.0000 0.0000 0.0000 33.0000 0.0000 18.4682 15.5318fval = 122.0000exitflag = 1 即运输方案为：甲市场的货由B厂送20箱；乙市场的货由A厂送35箱；丙商场的货由C厂送33箱；丁市场的货由B厂送18箱，再由C厂送16箱。
最低总运费为：122元 2 0-1整数规划求极值形如： min f T X s.t A X≤b Aeq X =beqxi为0或1 其中X为n维未知向量，f T=[f1,f2,…fn]为目标函数系数向量，小于等于约束系数矩阵A为m×n矩阵，b为其右端m维列向量，Aeq为等式约束系数矩阵，beq为等式约束右端常数列向量lb,ub为自变量取值上界与下界约束的n维常数向量2.1分支定界法在Matlab中提供了bintprog函数实现0-1型线性规划，采用的是分支定界法原理，其调用格式如下： x=bintprog(f,A,b)x=bintprog(f,A,b,Aeq,beq)[x,fval]=bintprog(…)[x, fval, exitflag]=bintprog(…)[x, fval, exitflag, output]=bintprog(…)[x, fval, exitflag, output, lambda]=bintprog(…)说明：x=bintprog(f,A,b)返回值x为最优解向量x=bintprog(f,A,b,Aeq,beq) 作有等式约束的问题若没有不等式约束，则令A=[ ]、b=[ ] 。
[x,fval]=bintprog(…) 左端 fval 返回解x处的目标函数值例6】求解以下问题：max z=x1+1.2x2+0.8x3st2.1x1+2x2+1.3x3<=50.8x1+x2<=5x1+2.5x2+2x3<=82x2<=8x1,x2,x3为0或1解答:首先，将其改变成0-1整数规划函数bintprog要求的标准形式：st2.1x1+2x2+1.3x3<=50.8x1+x2<=5x1+2.5x2+2x3<=82x2<=8x1,x2,x3为0或1Matlab求解：c=[-1,-1.2,-0.8];A=[2.1,2,1.3;0.8,1,0;1,2.5,2;0,2,0];b=[5;5;8;8];[x,fval]=bintprog(c,A,b);xfmax=-fval运行结果可得：x= 1 1 0fmax = 2.2000【例7】某快餐连锁经营公司有7个地点（A1，A2，---，A7）可以设立快餐店，由于地理位置因素，设立快餐店时必须满足以下要求：A1，A2，A3三个地点最多可选两个，A1和A5至少选取一个，A6和A7至少选取一个已知各个地点设立快餐店的投入和预计收益如表所示。
已知目前公司有650万元可以投资问怎样投资才能使公司预计收益最高？地点A1A2A3A4A5A6A7利润/万元101181215125投资/万元1031409515019316080解：这是一个选址问题首先引入0-1变量xi;xi=1,选择Ai地址；xi=0,不选择Ai地址则该问题的数学模型可以表示如下：max z=10x1+11x2+8x3+12x4+15x5+12x6+5x7st(1)103x1+140x2+95x3+150x4+193x5+160x6+80x7<=650(2)x1+x2+x3<=2(3)x4+x5>=1(4)x6+x7>=1(5)xi=0,1matlab求解程序如下：c=[-10 -11 -8 -12 -15 -12 -5];A=[103 140 95 150 193 160 80;1 1 1 0 0 0 0; 0 0 0 -1 -1 0 0;0 0 0 0 0 -1 -1];b=[650;2;-1;-1];[x,fval]=bintprog(c,A,b);xf=fval*(-1)运行程序得：x=[1;1;0;1;0;1;1],fval=502.2 枚举法除了采用分支定界法原理计算0-1整数规划外，还可以采用枚举法枚举法程序bintLp_E.mfunction [x,f]=bintLp_E(c,A,b,N)%[x,f]=bintLp_E(c,A,b,N):用枚举法求解下列0-1线性规划%min f=c'*x,s.t.A*x<=b,x的分量全为0或1%其中N表示约束条件A*x<=b中的前N个是等式，N=0时可以省略%返回结果x是最优解，f是最优解处的函数值if nargin<4 N=0;endc=c(:);b=b(:);[m,n]=size(A);x=[];f=abs(c')*ones(n,1);i=1;while i<=2^n B=de2bi(i-1,n)'; t=A*B-b; t11=find(t(1:N,:)~=0); t12=find(t(N+1:m,:)>0); t1=[t11;t12]; if isempty(t1) f=min([f,c'*B]); if c'*B==f x=B; end end i=i+1;end注意：以上程序需保存至搜索路径之下才可以调用。
例8】 求解下列0-1型整数线性规划解：分别采用二种算法计算如下：c=[3 -2 5];a=[1 2 -1;1 4 -1 ;1 1 0;0 4 1];b=[2;4;3;6];[x,fval]=bintprog(c,a,b)[xx,ffval]=bintLp_E(c,a,b)计算结果如下：x = 0 1 0fval = -2xx = 0 1 0ffval = -2【例9】A1加工B1,A2加工B3,A3加工B2配件时所需总时间最少,只需10小时.练习2:某公司有A1,A2,A3三项业务需要B1,B2,B3三位业务员处理,每个业务员处理业务的费用如表所示,其中业务员B2不能处理业务A1,问应指派何人去完成何项业务,使所需总费用最少?B1B2B3A11500不能处理800A21200900750A3900800900解:设xij表示第i项业务被第j位业务员处理，其中不能处理时可以认为费用非常高，比如999999元，则依题意可得如下模型：min z=1500x11+999999x12+800x13+1200x21+900x22+750x23+900x31+800x32+900x33stx11+x12+x13=1x21+x22+x23=1x31+x32+x33=1x11+x21+x31=1x12+x22+x32=1x13+x23+x33=1xij=0,1编写matlab程序：c=[1500 999999 800 1200 900 750 900 800 900];Aeq=[1 1 1 0 0 0 0 0 0; 0 0 0 1 1 1 0 0 0; 0 0 0 0 0 0 1 1 1; 1 0 0 1 0 0 1 0 0; 0 1 0 0 1 0 0 1 0; 0 0 1 0 0 1 0 0 1];beq=ones(6,1);[x,fval]=bintprog(c,[],[],Aeq,beq)[xx,ffval]=bintLp_E(c,Aeq,beq,6)运行程序可得：x=xx=[0 0 1 0 1 0 1 0 0];fval=ffval2600答案:A1被B3处理，A2被B2处理，A3被B1处理时总费用最少,只需2600元.3 整数规划求极值MATLAB关于整数规划没有内带的函数，目前关于整数规划的自编程序已编好，以下是其中之一：function[x,val,status]=ip(f,A,b,Aeq,beq,lb,ub,M,e)options=optimset('display','off'); bound=inf;[x0,val0]=linprog(f,A,b,Aeq,beq,lb,ub,[],options);[x,val,status,b]=rec(f,A,b,Aeq,beq,lb,ub,x0,val0,M,e,bound); function[xx,val,status,bb]=rec(f,A,b,Aeq,beq,lb,ub,x,v,M,e,bound)options=optimset('display','off');[x0,val0,status0]=linprog(f,A,b,Aeq,beq,lb,ub,[],options);if status0<=0|val0>bound xx=x;val=v;status=status0;bb=bound; return;endind=find(abs(x0(M)-round(x0(M)))>e);if isempty(ind) status=1; if val0 一般e=5.96e-08;该函数返回变量如下：X:整数规划的解Val;目标函数最优值Status=1 如果成功 =0 如果迭代到线性规划的最大迭代次数 =-1 如果没有解【例10】求解满足以下条件的极值：max z=20x1+10x2st5x1+4x2<=242x1+5x2<=13x1,x2>=0x1,x2取整数分析：由于函数是求极小值，故先变换成求负的最小值，然后变成正的最大值具体如下：c=[-20,-10];A=[5 4;2 5];b=[24;13];lb=[0 0];M=[1 2];e=5.96e-08;[x,val,status]=ip(c,A,b,[],[],lb,[],M,e);x,maxf=-val,status计算结果如下：x = 4 1maxf = 90.0000status = 1【例11】现有一个容积为36立方米,最大载重40吨的集装箱,需装入两种产品.产品甲为箱式包装,每箱体积0.3立方米,重0.7吨,每箱价值1.5万元;产品乙为袋式包装,每袋体积为0.5立方米,重0.2吨,每袋价值1万元.问箱子应当装入多少产品甲(不可拆开包装)以及多少产品乙(可以拆开包装)才能使集装箱载货价值最大?模型求解：假设应装入x1箱甲产品，x2袋乙产品，则问题变成如下数学模型：max z=1.5x1+x2st0.3x1+0.5x2<=360.7x1+0.2x2<=40x1,x2>=0x1取整数 编程如下：c=[-1.5,-1];A=[0.3 0.5;0.7 0.2];b=[36;40];lb=[0 0];M=[1];e=5.96e-08;[x,val,status]=ip(c,A,b,[],[],lb,[],M,e);x,fmax=-val,status计算结果为：x = 44.0000 45.6000fmax = 111.6000status = 1。
