Exp08姜启源数学建模课件第八章约束优化[高等教学]

1、大学数学实验大学数学实验 Mathematical Experiments 实验实验8 8 约束优化约束优化 1 优化问题三要素：优化问题三要素：决策变量决策变量；目标函数目标函数；约束条件约束条件 约约 束束 条条 件件决策变量决策变量 优化问题的一般形式优化问题的一般形式 n j i Dx ljxg mixhts xf ,.,1, 0)( ,.,1, 0)(. . )(min 当最优解在可行域边界上取得时当最优解在可行域边界上取得时 不能用无约束优化方法求解不能用无约束优化方法求解 目标函数目标函数 2 约束优化的分类约束优化的分类 线性规划线性规划(LP) 目标和约束均为线性函数目标和约

2、束均为线性函数 非线性规划非线性规划(NLP) 目标或约束中存在非线性函数目标或约束中存在非线性函数 二次规划二次规划(QP) 目标为二次函数、约束为线性目标为二次函数、约束为线性 整数规划整数规划(IP) 决策变量决策变量(部分部分)为整数为整数 整数整数线性线性规划规划(ILP) 整数整数非线性非线性规划规划(INLP) 0-1规划规划整数决策变量只取或整数决策变量只取或 n j i Dx ljxg mixhts xf ,.,1, 0)( ,.,1, 0)(. . )(min 连连 续续 优优 化化 离离 散散 优优 化化 数学规划数学规划 (Math. Prog.) 3 本实验基本内容本

3、实验基本内容 2. 基本原理和算法基本原理和算法 3. MATLAB实现实现 1. 问题与模型问题与模型 NLPNLP LPLP QPQP 连续优化连续优化 4 1桶桶 牛奶牛奶 3千克千克A1 12小时小时 8小时小时 4千克千克A2 或或 获利获利12元元/千克千克 获利获利8元元/千克千克 0.8千克千克B1 2小时小时,1.5元元 1千克千克 获利获利22元元/千克千克 0.75千克千克B2 2小时小时,1.5元元 1千克千克 获利获利16元元/千克千克 制订生产计划，使每天净利润最大制订生产计划，使每天净利润最大 15元可增加元可增加1桶牛奶，应否投资？桶牛奶，应否投资？ 50桶牛奶

4、桶牛奶, 480小时小时 至多至多100公斤公斤A1 B1，B2的获利经常有的获利经常有10%的波动，对计划有无影响？的波动，对计划有无影响？ 实例实例1: 奶制品生产销售计划奶制品生产销售计划 聘用临时工人增加劳动时间，工资最多每小时几元？聘用临时工人增加劳动时间，工资最多每小时几元？ 5 1桶桶 牛奶牛奶 3千克千克 A1 12小时小时 8小时小时 4千克千克 A2 或或 获利获利12元元/千克千克 获利获利8元元/千克千克 0.8千克千克 B1 2小时小时,1.5元元 1千克千克 获利获利22元元/千克千克 0.75千克千克 B2 2小时小时,1.5元元 1千克千克 获利获利16元元/千

5、克千克 出售出售x1 千克千克 A1, x2 千克千克 A2， x3千克千克 B1, x4千克千克 B2 原料原料 供应供应 劳动劳动 时间时间 加工能力加工能力 决策决策 变量变量 目标目标 函数函数 利润利润 约束约束 条件条件 非负约束非负约束 0, 61 xx x5千克千克 A1加工加工B1， x6千克千克 A2加工加工B2 654321 5 . 15 . 11622812xxxxxxzMax 50 43 6251 xxxx 48022 )(2)(4 65 6251 xx xxxx 100 51 xx 附加约束附加约束 53 80 x.x 64 750 x.x LP 6 50万元基金用

6、于投资三种股票万元基金用于投资三种股票A、B、C： A每股年期望收益每股年期望收益5元元(标准差标准差2元元)，目前市价，目前市价20元；元； B每股年期望收益每股年期望收益8元元(标准差标准差6元元)，目前市价，目前市价25元；元； C每股年期望收益每股年期望收益10元元(标准差标准差10元元)，目前市价，目前市价30元；元； 股票股票A、B收益的相关系数为收益的相关系数为5/24; 股票股票A、C收益的相关系数为收益的相关系数为0.5; 股票股票B、C收益的相关系数为收益的相关系数为0.25。 实例实例2：投资组合问题：投资组合问题 如期望今年得到至少如期望今年得到至少20%20%的投资回

7、报，应如何投资？的投资回报，应如何投资？ 投资回报率与风险的关系如何？投资回报率与风险的关系如何？ 假设：假设：1、基金不一定要用完（不用不计利息或贬值）、基金不一定要用完（不用不计利息或贬值） 2、风险通常用收益的方差或标准差衡量、风险通常用收益的方差或标准差衡量 7 决策变量决策变量 x1 、x2和和 x3 分别表示投资分别表示投资A、B、C的数量的数量 （国内股票通常以（国内股票通常以“一手一手”（100股）为最小单位出售，股）为最小单位出售， 这里以这里以100股为单位，期望收益以百元为单位）股为单位，期望收益以百元为单位） 实例实例2：投资组合问题：投资组合问题 A、B、C每手每手(

8、百股百股)的收益分别记为的收益分别记为S1,S2和和S3(百元百元)： ES1=5, ES2=8, ES3=10， DS1=4, DS2=36, DS3=100， r12=5/24, r13=-0.5，r23=-0.25 总收益总收益 S=x1S1+x2S2+x3S3 ：是一个随机变量：是一个随机变量 15),cov( 10),cov( 25),cov( 322332 311331 211221 DSDSrSS DSDSrSS DSDSrSS 8 投资风险（总收益的方差）为投资风险（总收益的方差）为 实例实例2：投资组合问题：投资组合问题 总期望收益为总期望收益为 Z1=ES= x1ES1+x

9、2ES2+x3ES3=5x1+8x2+10 x3 总收益总收益 S=x1S1+x2S2+x3S3 ：是一个随机变量：是一个随机变量 323121 2 3 2 2 2 1 32323131 21213 2 32 2 21 2 1 332233112211 3322113322112 30205100364 ),cov(2),cov(2 ),cov(2 ),cov(2),cov(2),cov(2 )()()()( xxxxxxxxx SSxxSSxx SSxxDSxDSxDSx SxSxSxSxSxSx SxDSxDSxDSxSxSxDZ 9 二次规划（二次规划（QP）模型）模型 实例实例2：投资

10、组合问题：投资组合问题 s.t. 5x1 +8x2+10 x3 1000 20 x1+25x2+30 x3 5000 x1，x2，x3 0 323121 2 3 2 2 2 12 30205100364minxxxxxxxxxZ 通过试探发现通过试探发现 从从0.00010.1以以0.0001的步长变化的步长变化 就可以得到很好的近似结果就可以得到很好的近似结果 Min Z =Z2 - Z1 s.t. 20 x1+25x2+30 x3 5000 x1，x2，x3 0 加权加权 模型模型 10 实例实例3 3：供应与选址：供应与选址 某公司有6个建筑工地，位置坐标为(ai, bi) (单位：公里

11、), 水泥日用量di (单位：吨） i a1.258.750.55.7537.25 b1.250.754.7556.57.75 d3547611 1)现有 2 料场，位于 A (5, 1), B (2, 7), 记(xj,yj),j=1,2, 日储量 ej各有 20 吨。 假设：假设：料场和工地之间有直线道路料场和工地之间有直线道路 目标：制定每天的供应计划，即从 A, B 两料场分别 向各工地运送多少吨水泥，使总的吨公里数最小。 11 2 , 1, 6,.,1, 0 2 , 1, 6,.,1,. . )()(min 6 1 2 1 2 1 6 1 2/122 jic jec idcts by

12、axc ij jij i iij j ji ijijij 线性规划模型线性规划模型 决策变量：决策变量： ci j 12维维 (料场料场j到到 工地工地i运量）运量） 12 2）改建两个新料场，需要确定新料场位置(xj,yj)和 运量cij ，在其它条件不变下使总吨公里数最小。 . 2 , 1, 6,.,1, 0 , 2 , 1, , 6,.,1,. . )()(min 6 1 2 1 2 1 6 1 2/122 jic jec idcts byaxc ij jij i iij j ji ijijij 决策变量：决策变量： ci j， ，(xj,yj)16维 维 非线性规划模型非线性规划模型

13、13 求解线性规划求解线性规划(LP)的基本原理的基本原理 基本模型基本模型 0,. . ,min)max( xbAxts Rxxczor nT mnmn RbRARc , 二维线性规划的图解法二维线性规划的图解法 一般线性规划的单纯形算法一般线性规划的单纯形算法 一般线性规划的内点算法一般线性规划的内点算法 14 二维线性规划的图解法二维线性规划的图解法 x2 x10 L1 L2 L3 z1=0 z2=2 z3=6 z4=10 z5=13 grad z x* * 起作用约束：起作用约束：L2;L3 最优解（最优解（4，1），最优值），最优值zmax=13 L1 L2 L3 2 21 xx 1

14、5 求解求解LP的基本思想的基本思想 从可行域的某一顶点开始，只需在有限多个顶点中从可行域的某一顶点开始，只需在有限多个顶点中 一个一个找下去，一定能得到一个一个找下去，一定能得到最优解最优解。 LP的约束和目标函数均为线性函数的约束和目标函数均为线性函数 2维维 可行域可行域 线段组成的凸多边形线段组成的凸多边形 目标函数目标函数 等值线为直线等值线为直线 最优解最优解 凸多边形的某个顶点凸多边形的某个顶点 n维维 超平面组成的凸多面体超平面组成的凸多面体 等值线是超平面等值线是超平面 凸多面体的某个顶点凸多面体的某个顶点 算法：怎样从一点转到下一点，尽快找到最优解。算法：怎样从一点转到下一

15、点，尽快找到最优解。 16 x2 x10 L1 L2 L3 x2 x10 L1 L2 L3 x2 x10 L1 L2 求解求解LPLP的特殊情形的特殊情形 L1 L2 L3 无可行解无可行解无最优解无最优解(无界无界)最优解不唯一最优解不唯一 321 4123Lxx 3 L无 x2 x10 L1 L2 L3 z=c 2 21 xx 321 413Lxx 17 线性规划的标准形和基本性质线性规划的标准形和基本性质 nmRA xbAxts xcz nm T , 0,. min 标标 准准 形形 加入松弛变量加入松弛变量/剩余变量将不等式变为等式剩余变量将不等式变为等式 0, 0,1423 0, 2

16、2 0003min 21 5521 4421 54321 xx xxxx xxxx xxxxxz 0, 2. 3321 xxxxts 2 21 xx 一般还假定一般还假定 b0 rank(A)=m 18 0, ,1423 , 22 2. . 0003min 54321 521 421 321 54321 xxxxx xxx xxx xxxts xxxxxz 0, ,1423 , 22 2 54321 521 421 321 xxxxx xxx xxx xxx 对标准形求解对标准形求解 先求可行解先求可行解 nmRA xbAxts xcz nm T , 0,. min 0,xbAx 满足 再在有

17、限个可行解中寻找最优解再在有限个可行解中寻找最优解 19 1423 22 2 521 421 321 xxx xxx xxx 10023 01021 00111 A 可逆， BNBA, T NB T xxx, bNxBxAx NB 54321 ppppp bBxx BN 1 0 x2 x1O T xpppB)14, 2 , 2 , 0 , 0( 543 T xpppB) 8 , 0 , 4 , 0 , 2( 531 T xpppB)16, 0 , 3 , 1, 0( 532 O点点 Q点点 R点点 Q R 基本可行解基本可行解x ：Ax=b, x 0 (xB 0，xN=0) T xpppB)0

18、 , 0 , 5 , 1 , 4( 321 P点点 P B: 基基(矩阵矩阵) x: 基基(本本)解解 20 可行域存在时，必是凸多面体可行域存在时，必是凸多面体(可能无界可能无界)； 最优解存在时，必在可行域的顶点取得最优解存在时，必在可行域的顶点取得(可能无界可能无界)； 基本可行解对应于可行域的顶点基本可行解对应于可行域的顶点(极点极点)。 LPLP基本性质基本性质 最优解只需在有限个可行解（基本可行解）中寻找最优解只需在有限个可行解（基本可行解）中寻找 单纯形法的基本思路单纯形法的基本思路 从一个顶点转换到另一个顶点（即构成从一个顶点转换到另一个顶点（即构成xB和和xN 的向量的向量p

19、i间的转换），使目标函数逐步下降。间的转换），使目标函数逐步下降。 LP的通常解法是单纯形法的通常解法是单纯形法(G. B. Dantzig, 1947) 21 内点算法内点算法(Interior point method) 20世纪世纪80年代人们提出的一类新的算法年代人们提出的一类新的算法内点算法内点算法 也是迭代法，但不再从可行域的一个顶点转换到另一个也是迭代法，但不再从可行域的一个顶点转换到另一个 顶点，而是直接从可行域的内部逼近最优解。顶点，而是直接从可行域的内部逼近最优解。 LPLP其他算法其他算法 有效集有效集(Active Set)方法方法 LP是是QP的特例（只需令所有二次项

20、为零即可）的特例（只需令所有二次项为零即可） 可以用可以用QP的算法解的算法解LP(如如: 有效集方法，详见下节有效集方法，详见下节) 22 x,fval,exitflag,output,lambda = linprog(c,A1,b1,A2,b2,v1,v2,x0,opt) MATLAB MATLAB 求解求解 LPLP 212211 ,. . min vxvbxAbxAts xcz T 输入：输入： x0 x0初始解初始解( (缺省时一般为缺省时一般为0)0) opt opt MATLAB MATLAB控制参数控制参数 中间所缺参数项补中间所缺参数项补 输出：输出： lambda lamb

21、da LagrangeLagrange乘子，乘子，维维 数等于约束个数，非零分量对数等于约束个数，非零分量对 应于起作用约束应于起作用约束 lambda.ineqlin: 对应 A1xb1 lambda. eqlin : 对应 A2x = b2 lambda. lower : 对应 v1 x lambda. upper : 对应 x v2Exam0802.m 23 MATLAB MATLAB 求解求解 LPLP opt MATLAB控制参数:三种三种算法选择 缺省时采用大规模算法（是一种内点算法）；缺省时采用大规模算法（是一种内点算法）； 当当opt中中“LargeScale”参数设置为参数设

22、置为“off”时，采用中规模算法，时，采用中规模算法， 该模式下缺省的算法是二次规划的算法（一种有效集方法）；该模式下缺省的算法是二次规划的算法（一种有效集方法）； 当当opt中中“LargeScale”参数设置为参数设置为“off”，并且，并且“Simplex”参参 数设置为数设置为“on”时，采用单纯形算法。时，采用单纯形算法。 只有有效集方法可以由用户提供初始解只有有效集方法可以由用户提供初始解x0，其他两个算法则不，其他两个算法则不 需要（即使提供了也会被需要（即使提供了也会被MATLAB忽略）。忽略）。 Exam0801.m x,fval,exitflag,output,lambda

23、 = linprog(c,A1,b1,A2,b2,v1,v2,x0,opt) 24 c=12 8 22 16 -1.5 -1.5; A1=4 3 0 0 4 3;2 1 0 0 3 2;1 0 0 0 1 0; b1=600 240 100; A2=0 0 1 0 -0.8 0;0 0 0 1 0 -0.75; b2=0 0; v1=0 0 0 0 0 0; x,z0,ef,out,lag=linprog(-c,A1,b1,A2,b2,v1) lag.ineqlin, lag.eqlin 实例实例1: 奶制品生产销售计划奶制品生产销售计划 654321 5 . 15 . 11622812xxx

24、xxxzMax 50 43 6251 xxxx 48022)(2)(4 656251 xxxxxx 100 51 xx 08 . 0 53 xx 075. 0 64 xx 0 654321 x ,x ,x ,x ,x ,x 6003434 6521 xxxx 240232 6521 xxxx x=(0,168,19.2,0,24,0) ; z = -z0 =1730.4; lag.ineqlin =(1.58;3.26; 0.00) ; 25 15元可增加元可增加1桶牛奶，应否投资？桶牛奶，应否投资？ 实例实例1: 奶制品生产销售计划奶制品生产销售计划 654321 5 . 15 . 1162

25、2812xxxxxxzMax 50 43 6251 xxxx 48022)(2)(4 656251 xxxxxx 100 51 xx 53 80 x.x 64 750 x.x 0 654321 x ,x ,x ,x ,x ,x x=(0,168,19.2,0,24,0) ; z = -z0 =1730.4 lag.ineqlin =(1.58;3.26; 0.00) ; 6003434 6521 xxxx 601 z1=1731.98 dz=z1-z = 1731.98-1730.4=1.58 dz=lag.ineqlin(1) dz*12=1.58*12= 18.9615 应该投资！应该投资

26、！ “影子价格影子价格” 26 实例实例1: 奶制品生产销售计划奶制品生产销售计划 聘用临时工人增加劳动时间，工资最多每小时几元？聘用临时工人增加劳动时间，工资最多每小时几元？ 654321 5 . 15 . 11622812xxxxxxzMax 50 43 6251 xxxx 48022)(2)(4 656251 xxxxxx 100 51 xx 53 80 x.x 64 750 x.x 0 654321 x ,x ,x ,x ,x ,x x=(0,168,19.2,0,24,0) ; z = -z0 =1730.4 lag.ineqlin =(1.58;3.26; 0.00) ; 6003

27、434 6521 xxxx 240232 6521 xxxx lag.ineqlin(2)=3.26， 所以1小时劳动时间的影 子价格应为3.26/2=1.63， 即单位劳动时间增加的利 润是1.63(元) 27 B1，B2的获利经常有的获利经常有10%的波动，对计划有无影响？的波动，对计划有无影响？ 实例实例1: 奶制品生产销售计划奶制品生产销售计划 654321 5 . 15 . 11622812xxxxxxzMax 50 43 6251 xxxx 48022)(2)(4 656251 xxxxxx 100 51 xx 53 80 x.x 64 750 x.x 0 654321 x ,x

28、,x ,x ,x ,x x=(0,168,19.2,0,24,0) ; z = -z0 =1730.4 lag.ineqlin =(1.58;3.26; 0.00) ; 若每公斤若每公斤B1的获利下降的获利下降10%, 应将目标函数中应将目标函数中x3的系数改的系数改 为为19.8，重新计算发现最优，重新计算发现最优 解和最优值均发生了变化解和最优值均发生了变化 若若B2的获利向上波动的获利向上波动10%， 原计划也不再是最优的原计划也不再是最优的 MATLAB没有给出这种没有给出这种 敏感性分析敏感性分析的结果的结果 (LINDO/LINGO可以可以) 28 非线性规划非线性规划(NLP)(

29、NLP)基本原理基本原理 ljxgts xfz j x n ,.,1, 0)(. . )(min g1 d g2=0 g1=0 g3=0 gj0 o x 设设x为可行解，为可行解， 位于约束边界位于约束边界 1 , 0)(Jjxg j J1起作用约束起作用约束(j=1) 2 , 0)(Jjxg j J2不起作用约束不起作用约束(j=2,3) )0 ( 0 Gdx 可行方向可行方向 下降方向下降方向 )0()()( 0 xfdxf (G) )1 ( )( 0)( 1 Jj dxg T j )2(0)(dxf T )()()()( 2 Odxgxgdxg T jjj 不等式约束不等式约束 29 x

30、为最优解为最优解 不存在满足不存在满足(1),(2)的的d ) 1 ()(0)( 1 Jjdxg T j )2(0)(dxf T g3(x)=0 g1(x)=0 g2(x)=0 G(gj0) x O g1(x) g2(x) f(x) x0 g1(x0) d f(x)=c c g1(x0) 若若x沿沿d方向既可行又下降，则方向既可行又下降，则x不是最优解不是最优解 观察观察 f 的梯度能表示成的梯度能表示成g1 和和g2的梯度的非负的梯度的非负 线性组合线性组合! 最优解的必要条件最优解的必要条件 30 0)()( 1 xgxf j l j i ljxg jj , 2 , 1, 0)( KKT条

31、件条件 互补性条件 ljxg mixhts xf j i x n ,.,1,0)( ,.,1,0)(. )(min ，和则存在 线性无关， 0 )(, 1 ji ji Jjgh 0)()()( 11 xgxhxf j l j i m i ii ljxg jj , 1, 0)( 0, 1 l )()( 1 Jjxg j 且且线性无关，则存在线性无关，则存在 若若x为最优解，为最优解， 最优解的必要条件最优解的必要条件 31 KKTKKT条件的几何解释条件的几何解释 Q P f - g1- g2 f - g1 - g2 0)( 04)( 010)(. . ) 3() 7()(min 23 212

32、2 2 2 11 2 2 2 1 xxg xxxg xxxgt s xxxf 0)()( 1 xgxf j l j i ljxg jj , 2 , 1, 0)( x2 x1 0 最优解在最优解在P(3,1)取得取得 0)(2)()( , 0, 2, 1 21 321 PgPgPf P(3,1)是是KKT点点 其它点其它点(如如Q)均不是均不是 (7,3) 32 二次规划二次规划(QP)(QP)及有效集方法及有效集方法 bAxts cxHxxxf T . . 2 1 )(min 当H为对称阵，称二次规划 当H正定时，称凸二次规划 凸二次规划性质： 最优点KKT点； 局部最优解全局最优解； mTT

33、 RbAxcxHxxxL),( 2 1 ),( 最优解方程最优解方程 0 0 bAx AcHx TT L函数函数 等式约束 下 的Lagrange乘子法 bAx 33 解二次规划的有效集方解二次规划的有效集方法法 基本思想：基本思想：对于不等式约束的二次规划，在某可行点 处将不起作用约束去掉，起作用约束视为等式约束起作用约束视为等式约束， 通过求解等式约束的二次规划来改进可行点。 若x为(1)的最优解，则它也是(2)的最优解 若x为(1)的可行解，又是(2)的KKT点，且 L乘子非负，则它必是(1)的最优解。 bAxts cxHxxxf T . . 2 1 )(min (1) 1 ,. . 2

34、 1 )(min Jjbxat s cxHxxxf jj T )(列的第是jAa j (2) 基本基本 原理原理 34 *, 0. . )*()*()*( 2 1 )*(min Jjdats dxcdxHdxdxf j T 若d*0，且x*+d*可行则继续;否则确定步长*( 1)使 ap(x*+*d*)=bp，pJ*，则有效集修正为J*p。 设(1)的可行点为x*，有效集 记作J*，用L乘子法求解: 基本步骤基本步骤 得d*, * 若d*=0, 则x*为(2)最优解; 当 *非负时x*是(1)最优解 若d*=0, 且( *)q 1 g = . % gradient of the functio

35、n if nargout 2 H = . % Hessian end 41 nlcon.m给出约束,GradConstr=on时还给出梯度,形式为 212211 21 , , 0)(, 0)(. . )(min vxvbxAbxA xcxcts xfz function c1,c2,GC1,GC2 = nlcon(x) c1 = . % nonlinear inequalities at x c2 = . % nonlinear equalities at x if nargout 2 GC1 = . % gradients of c1 GC2 = . % gradients of c2 en

36、d MATLABMATLAB求解求解 约束约束NLPNLP 例例4 Exam0804.m ) 12424(min 221 2 2 2 1 1 xxxxxe x 010, 05 . 1 212121 xxxxxx 01 2 2 1 xx 42 MATLABMATLAB优化工具箱优化工具箱能求解的优化模型能求解的优化模型 优化工具箱优化工具箱3.0 (MATLAB 7.0 R14) 连续优化连续优化离散优化离散优化 无约束优化无约束优化 非线性非线性 极小极小 fminunc 非光滑非光滑(不可不可 微微)优化优化 fminsearch 非线性非线性 方程方程(组组) fzero fsolve 全局全局 优化优化 暂缺暂缺 非线性非线性 最小二乘最小二乘 lsqnonlin lsqcurv