(大学数学实验 约束优化问题一般形式 大学数学实验 f(x) 当最优解在可行域边 界上取得时,不能用 Experiments in Mathematics 8(x)≤0,j=1,1无约束优化方法求解 实验8约束优化 约束优化基本内客:LP,QP,NLP 问题与模型 半大总数科系 2.基本原理 3.算法和 MATLAB实现 实例1:奶制品生产销售计划 2小时千一我利2元千克 牛奶或 口时元[千览一获利2元千克 0.8千克B1一获利2元千克 8小时奸千克A一获利元千克 2小时,1.5元 2小时15元1075千宽B→获利16元千克 8小时4千克→获利元千克 决算出售x1千克A1,x2千克A2,x千克B1,x千克B2 50桶牛奶480小时小解1元03千克B]一获利16元千克 变量x千克A加工B,x千克A加工B2 至多100公斤A1制订生产计划,使每天净利润最大 画数利润 + =12x1+8x2+22x2+16x4-1.5x5-1.5x6 15元可增加1桶牛奶,应否投资? 加工能力x1+x5≤1 约束 聘用临时工人增加劳动时间,工资最多每小时几元? 条件 4(x+x)+2(x2+)的取x3=0.8 x4=0.75x6 B1,B2的获利经常有10%的波动,对计划有无影响? +2x+2x6≤480非负约東x1,…x620 (学学奖 实例2:投资组合问题 实例2:投资组合问题 50万元基金用于投资三种股票A、B、C: A、B、C每手(百股)的收益分别记为S1S2和S百元 B每股年期望收益8元(标准差6元),目前市价25元; ES1=5,ES2=8,ES3=10, N(S1,S)=h√DS√DS2=25 年期望收益10元(标准差10元),目前市价30元 股票A、B收益的相关系数为5/4 DS=4,DS2=36,DS3=100, 股票A、C收益的相关系数为0.5; r12=5/24,r1=0.5,r2=0.25co(s2,S3)=53√DS2√DS3 股票B、C收益的相关系数为0.25。 决策变量x1、x2和x3分别表示投资A、B、C的数量 如期望今年得到至少20%的投资回报,应如何投资? (国内股票通常D 100股)为最小单位出售 投资回报率与风险的关系如何? 这里以100股为单位,期望收益以百元为单位 假设:1、基金不一定要用完(不用不计利息或贬值) 2、风险通常用收益的方差或标准差衡量 总收益S=x1S1+xS2+xS3:是一个随机变量
1 大学数学实验 Experiments in Mathematics 实验8 约束优化 清华大学数学科学系 约束优化基本内容:LP,QP,NLP 2. 基本原理 3. 算法和MATLAB实现 1. 问题与模型 约束优化问题一般形式 g x j l s t h x i m f x j i x n ( ) 0, 1,..., . . ( ) 0, 1,..., min ( ) ≤ = = = ∈ℜ 当最优解在可行域边 界上取得时，不能用 无约束优化方法求解 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桶牛奶, 480小时 至多100公斤A1 • B1，B2的获利经常有10%的波动，对计划有无影响？ 实例1: 奶制品生产销售计划 • 聘用临时工人增加劳动时间，工资最多每小时几元？ 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元/千克 出售x1 千克 A1, x2 千克 A2， x3千克 B1, x4千克 B2 原料 供应 劳动 时间 加工能力 决策 变量 目标 函数 利润 约束 条件 非负约束 , 0 x1 Lx6 ≥ x5千克 A1加工B1， x6千克 A2加工B2 1 2 3 4 5 6 Max z =12x +8x + 22x +16x −1.5x −1.5x 50 3 4 1 5 2 6 ≤ + + x + x x x 2 2 480 4( ) 2( ) 5 6 1 5 2 6 + + ≤ + + + x x x x x x 100 x1 + x5 ≤ 附加约束 3 5 x = 0.8x 4 6 x = 0.75x 50万元基金用于投资三种股票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%的投资回报，应如何投资？ •投资回报率与风险的关系如何？ 假设：1、基金不一定要用完（不用不计利息或贬值） 2、风险通常用收益的方差或标准差衡量 决策变量 x1 、x2和 x3 分别表示投资A、B、C的数量 （国内股票通常以“一手”（100股）为最小单位出售， 这里以100股为单位，期望收益以百元为单位） 实例2：投资组合问题 A、B、C每手(百股)的收益分别记为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 ：是一个随机变量 cov( , ) 15 cov( , ) 10 cov( , ) 25 2 3 23 2 3 1 3 13 1 3 1 2 12 1 2 = = − = = − = = S S r DS DS S S r DS DS S S r DS DS
实例2:投资组合问题 实例2:投资組合问题 总收益S=xS1+x2S2+x353:是一个随机变量 minZ2=4x2+36x2+100x2+5xx2-20xx3-30x2x3 总期望收益为 st.5x;+8x,+10x21000 ZIES=x ES,+x,ES,+x3ES,=5x,+&r,+10x 0x1+25x2+30x3≤5000 投资风险(总收益的方差)为 加权 Min Z= BZ,-Z Z2=D(,S+x2S2+x,S,)=D(r S1)+D( S2)+D(r,,) 模型 st.20x1+25x2+30x3≤5000 +2cov(r S, x2S2)+2cov(r S, x3 S3)+2 cov(r2S,,x3 S3 =x DS,+x2DS,+xfDS, +2x, cov(S,,S,) 通过试探发现B从0.00101以0.0001的步长变 化就可以得到很好的近似结果 +2x, cov(S,,S,)+2x2Y3 cov(S2,S,) =4x+36x2+1003+5xx2-20xx3-390为9e 二次规划(QP)模型 某公司有个建筑地,位置坐标为 实例3:供应与选址 决簟变量:c +(yj 水泥日用量d(单位:吨) (料场到工地 运量)-12维 0.754.75 l)现有2料场,位于A(5,1),B(2,7) 记(x,y)=1,2,日储量e各有20吨。 假设:料场和工地之间有直线道路 目标:制定每天的供应计划,即从A,B两料场分别 线性規划模型 向各工地运送多少吨水泥,使总的吨公里数最小 (学学奖 2)改建两个新料场,需要确定新料场位置(xy)和 约束优化问题的分类 运量,在其它条件不变下使总吨公里数最小 cyxy-a2+(y-b)22决策变量 线性规划P)目标和约束均为线性函数 c(xy)16维 ·非线性规划(NLP目标或约束中存在非线性函数 √二次规划(QP目标为二次函数、约束为线性 cy=d,i=1-6 整数规划(P)决策变量部分)为整数 非线性规划模型 √整数线性规划(LLP) ≤2 √整数非线性规划(INLP
2 投资风险（总收益的方差）为 实例2：投资组合问题 总期望收益为 Z1=ES= x1ES1+x2ES2+x3ES3=5x1+8x2+10x3 总收益 S=x1S1+x2S2+x3S3 ：是一个随机变量 1 2 1 3 2 3 2 3 2 2 2 1 1 3 1 3 2 3 2 3 3 1 2 1 2 2 2 3 2 1 2 2 1 1 1 2 2 1 1 3 3 2 2 3 3 2 1 1 2 2 3 3 1 1 2 2 3 3 4 36 100 5 20 30 2 cov( , ) 2 cov( , ) 2 cov( , ) 2cov( , ) 2cov( , ) 2cov( , ) ( ) ( ) ( ) ( ) x x x x x x x x x x x S S x x S S x DS x DS x DS x x S S x S x S x S x S x S x S Z D x S x S x S D x S D x S D x S = + + + − − + + = + + + + + + = + + = + + 二次规划（QP）模型 实例2：投资组合问题 s.t. 5x1 +8x2+10x3 ≥ 1000 20x1+25x2+30x3 ≤ 5000 x1，x2，x3 ≥ 0 1 2 1 3 2 3 2 3 2 2 2 2 1 min Z = 4x + 36x +100x + 5x x − 20x x − 30x x 通过试探发现 β从0.0001~0.1以0.0001的步长变 化就可以得到很好的近似结果 Min Z =βZ2 - Z1 s.t. 20x1+25x2+30x3 ≤ 5000 x1，x2，x3 ≥ 0 加权 模型 实例3：供应与选址 某公司有6个建筑工地，位置坐标为(ai , bi ) (单位：公里), 水泥日用量di (单位：吨） i １ ２ ３ ４５ ６ a 1.25 8.75 0.5 5.75 3 7.25 b 1.25 0.75 4.75 5 6.5 7.75 d 3 5 4 7 6 11 1)现有 2 料场，位于 A (5, 1), B (2, 7), 记(xj ,yj ),j=1,2, 日储量 ej各有 20 吨。 假设：料场和工地之间有直线道路 目标：制定每天的供应计划，即从 A, B 两料场分别 向各工地运送多少吨水泥，使总的吨公里数最小。 , 1,2 . . , 1,...,6 min [( ) ( ) ] 6 1 2 1 2 1 6 1 2 2 1/ 2 ≤ = = = − + − ∑ ∑ ∑∑ = = = = c e j s t c d i c x a y b ij j i ij i j j i ij j i j i 线性规划模型 决策变量：ci j (料场j到工地i的 运量）~12维 2）改建两个新料场，需要确定新料场位置(xj ,yj )和 运量cij ，在其它条件不变下使总吨公里数最小。 , 1,2 . . , 1,...,6 min [( ) ( ) ] 6 1 2 1 2 1 6 1 2 2 1/ 2 ≤ = = = − + − ∑ ∑ ∑∑ = = = = c e j st c d i c x a y b ij j i ij i j j i ij j i j i 决策变量： ci j，(xj ,yj )~16维 非线性规划模型 约束优化问题的分类 • 线性规划(LP) 目标和约束均为线性函数 • 非线性规划(NLP) 目标或约束中存在非线性函数 9 二次规划(QP) 目标为二次函数、约束为线性 • 整数规划(IP) 决策变量(部分)为整数 9 整数线性规划(ILP) 9 整数非线性规划(INLP)
求解线性规划(LP)的基本原理 二维线性规划的图解法 max(or minE= x,xER" max =3x +x s.1.x1-x2≥ Ax≤b.x≥0 rad z c∈R",A∈Rm,b∈R 3x,+2x:≤14~L 盘肌划的圆油 起作用约束:L2;L3 限划的单形 朝的内成油 最优解(4,1),最优值z=13 LP的约来和目标函敷均为线性函数 求解LP的特殊情形max=3x,+x x1-2x:≤2-L2 可行域线段组成的凸多边形超平面组成的凸多面体 3x,+2x.≤14-L3 目标函数等值线为直线 等值线是超平面 最优解凸多边形的某个顶点凸多面体的某个顶点 -3x1+2x214~L无L3 3x1+x2≤14~L3 求解LP的基本思想 从可行域的某一顶点开始,只需在有限多个顶点中 个一个找下去,一定能得到最优解 算法:怎样从一点转到下一点,尽快找到最优解。 无可行解无最优解(无界)最优解不唯 (学学奖 大学费学买验 线性规划的标准形和基本性质 对标准形求解mn=-3x-x2+0x2+0x+0x nz=-3x1-x2+0·x3+0x4+0 .x1-x2≥-2 mmn二=cx s..-x+x2+x3=2x≥0 x-2x2+x4=2, 3. ∈ 3x1+2x2+x5=14,x≥0 x1,x2≥0 先求可行解 x+x2+x3=2 加入松弛变量剩余变量将不等式变为等式 标 lIn 2 准s1.Ax=b,x≥0 形 A∈Rm×n 再在有限个可行解中寻找最优解
3 求解线性规划(LP)的基本原理 基本模型 . . , 0 max( min) , ≤ ≥ = ∈ st Ax b x or z c x x R T n n m n m c ∈ R A ∈ R b ∈ R × , , •二维线性规划的图解法 •一般线性规划的单纯形算法 •一般线性规划的内点算法 二维线性规划的图解法 x2 0 x1 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 , 0 3 2 14 2 2 . . 2 max 3 1 2 1 2 1 2 1 2 1 2 ≥ + ≤ − ≤ − + ≤ = + x x x x x x s t x x z x x x1 − x2 ≥ −2 求解LP的基本思想 从可行域的某一顶点开始，只需在有限多个顶点中 一个一个找下去，一定能得到最优解。 LP的约束和目标函数均为线性函数 2维 可行域 线段组成的凸多边形 目标函数 等值线为直线 最优解 凸多边形的某个顶点 n维 超平面组成的凸多面体 等值线是超平面 凸多面体的某个顶点 算法：怎样从一点转到下一点，尽快找到最优解。 x2 x1 0 L1 L2 L3 x2 x1 0 L1 L2 L3 x2 x1 0 L1 L2 求解LP的特殊情形 ~L1 ~L2 ~L3 无可行解 无最优解(无界) 最优解不唯一 1 2 ~ 3 − 3x + 2x ≥ 14 L 无 L3 x2 x1 0 L1 L2 L3 z=c , 0 3 2 14 2 2 . . 2 max 3 1 2 1 2 1 2 1 2 1 2 ≥ + ≤ − ≤ − + ≤ = + x x x x x x s t x x z x x 2 x1 − x2 ≥ − 3 1 2 14 ~ L3 x + x ≤ 线性规划的标准形和基本性质 A R m n s t Ax b x z c x m n T ∈ ≤ = ≥ = × , . . , 0 标 min 准 形 加入松弛变量/剩余变量将不等式变为等式 , 0 3 2 14, 0 2 2, 0 min 3 0 0 0 1 2 1 2 5 5 1 2 4 4 1 2 3 4 5 ≥ + + = ≥ − + = ≥ = − − + ⋅ + ⋅ + ⋅ x x x x x x x x x x z x x x x x s.t. −x1 +x2 +x3 =2, x3 ≥0 , 0 3 2 14 2 2 . . 2 max 3 1 2 1 2 1 2 1 2 1 2 ≥ + ≤ − ≤ − + ≤ = + x x x x x x s t x x z x x x1 − x2 ≥ −2 , , , , 0 3 2 14, 2 2, . . 2 min 3 0 0 0 1 2 3 4 5 1 2 5 1 2 4 1 2 3 1 2 3 4 5 ≥ + + = − + = − + + = = − − + ⋅ + ⋅ + ⋅ x x x x x x x x x x x st x x x z x x x x x , , , , 0 3 2 14, 2 2, 2 1 2 3 4 5 1 2 5 1 2 4 1 2 3 ≥ + + = − + = − + + = x x x x x x x x x x x x x x 对标准形求解 先求可行解 A R m n s t Ax b x z c x m n T ∈ ≤ = ≥ = × , . . , 0 min Ax = b, x ≥ 0 满足 再在有限个可行解中寻找最优解
A→[B,NB可逆 A=1-2010 LP基本性质 x7→[xB,x 3200 可行域存在时,必是凸多面体 3n+22+5=14=[BB2B2P2B Ax= BxB+ NxN=b 可行解对应于可行域的顶点 最优解存在时,必在可行域的顶点取得 B=IP3P4ps]→x=(00,2,214)7O点 最优解只需在有限个可行解中寻找 B=[乃351→x=(20,4,08)Q点 LP的遢常解油是单施形油 B=[P2P3P5]→x=(0-13016)R点 单形法的本尾隐是 B=[BP2P]→x=(41.50)P点 从一个顶点转换到另一个顶点(即构成xg和x的 基本可行解x:Ac=b,x20(x20,x=00⊙ 向量p间的转换),使目标函数下降最多 LP其他算法 MATLAB求解LP min2三cx 内点法( (Interior point method) S!.Ax≤b1,A2x=b2,V1≤x≤v 0世纪80年代人们提出的一类新的算法—内点算 [x, fval, exitflag, output, lambda] 法也是迭代法,但不再从可行域的一个顶点转换到另 linprog(c, Al, bl, A2, b2, vl, v2, x0, options) 个顶点,而是直接从可行域的内部逼近最优解 1 ambda~ lagrange乘子 有效集( Active Set)方法 x0~初始解(缺省时为0) 维数等于约束个数,非零分 options~ MATLAB控制参数 量对应于起作用的束 线性规划是二次规划的一个特例(只需令所有二次 项为零即可),因此也可以用二次规划的算法(如有效 中间所缺参数项补[ 集方法,详见下节关于二次规划的介绍)解线性规划 lambda upper lambda. lower (学学奖 大学酸学实验) MATLAB求解LP 实例1:奶制品生产销售计划 options~ MATLAB控制参数:三种算法选择 Max==12x1+8x2+22x3+16x4-1.5x5-1.5x 缺省时采用大规模算法〔是一种内点算法 曰4x1+3r2 +4x;+3x6 e”参数设为“om时,采用中规模算法 4x+x,)+2x2+x)+2x+2x.5480中2x+x +3xs+2x≤ 该模式下缺省的算法是二次规划的算法(一种有效集方法); 叫1282216-1.5-1.5 当opr中“ Largescale"参数设置为“om”,并且“ Simplex”参数设 l=430043;210032;100010; 量为“on"时,采用单纯形算法 x3-0.8xs=0 bl=60024010 A2=|001040.80:0001040.75l 只有有效集方法可以由用户提供初始解x0,其他两个算法则不 x,x2,x,x,x,x20b2=100 需要(即使提供了也会被忽略 vl=000000 X, zO, ef, out, lag linprog(-c, Al, bl, A2, b2, vl) ag.ineqlin, lag eqi Exam080l m =(0,168,192,0,24,0);二=-z0=1730 ag. ineqlin=(1.58:3.26;0.00
4 3 2 14 2 2 2 1 2 5 1 2 4 1 2 3 + + = − + = − + + = x x x x x x x x x ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ − − = 3 2 0 0 1 1 2 0 1 0 1 1 1 0 0 A A⇒[B,N]，B可逆 [ , ] T N T B T x ⇒ x x Ax = BxB + NxN = b [ ] = p1 p2 p3 p4 p5 xN xB B b 1 0 − = ⇒ = x2 x O 1 T B [ p p p ] x (0,0,2,2,14) = 3 4 5 ⇒ = T B [ p p p ] x (2,0,4,0,8) = 1 3 5 ⇒ = T B [ p p p ] x (0, 1,3,0,16) = 2 3 5 ⇒ = − O点 Q点 R点 Q R 基本可行解x ：Ax=b, x ≥0 (xB≥0，xN=0) T B [ p p p ] x (4,1,5,0,0) = 1 2 3 ⇒ = P点 P 可行域存在时，必是凸多面体； 可行解对应于可行域的顶点； 最优解存在时，必在可行域的顶点取得。 LP基本性质 最优解只需在有限个可行解中寻找 单纯形法的基本思路是 从一个顶点转换到另一个顶点（即构成xB和xN的 向量pi 间的转换），使目标函数下降最多。 LP的通常解法是单纯形法 内点法(Interior point method) 20世纪80年代人们提出的一类新的算法——内点算 法也是迭代法，但不再从可行域的一个顶点转换到另一 个顶点，而是直接从可行域的内部逼近最优解。 LP其他算法 有效集(Active Set)方法 线性规划是二次规划的一个特例（只需令所有二次 项为零即可），因此也可以用二次规划的算法（如有效 集方法，详见下节关于二次规划的介绍）解线性规划 [x,fval,exitflag,output,lambda] = linprog(c,A1,b1,A2,b2,v1,v2,x0,options) MATLAB 求解 LP 1 1 2 2 1 2 . . , , min s t A x b A x b v x v z c xT ≤ = ≤ ≤ = 输入： x0~初始解(缺省时为0) options ~ MATLAB控制参数 中间所缺参数项补[] 输出： lambda ~ Lagrange乘子， 维数等于约束个数，非零分 量对应于起作用约束 • lambda.ineqlin • lambda. eqlin • lambda. upper • lambda. lower Exam0802.m MATLAB 求解 LP options ~ MATLAB控制参数:三种算法选择 • 缺省时采用大规模算法（是一种内点算法）； • 当opt中“LargeScale”参数设置为“off”时，采用中规模算法， 该模式下缺省的算法是二次规划的算法（一种有效集方法）； • 当opt中“LargeScale”参数设置为“off”，并且“Simplex”参数设 置为“on”时，采用单纯形算法。 只有有效集方法可以由用户提供初始解x0，其他两个算法则不 需要（即使提供了也会被忽略）。 Exam0801.m 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: 奶制品生产销售计划 1 2 3 4 5 6 Max z = 12x + 8x + 22x +16x −1.5x −1.5x 50 3 4 1 5 2 6 ≤ + + x + x x x 4(x1 + x5 ) + 2(x2 + x6 ) + 2x5 + 2x6 ≤ 480 x1 + x5 ≤ 100 0.8 0 x3 − x5 = x4 − 0.75x6 = 0 0 x1 ,x2 ,x3 ,x4 ,x5 ,x6 ≥ 4x1 + 3x2 + 4x5 + 3x6 ≤ 600 2 3 2 240 x1 + x2 + x5 + x6 ≤ x=(0,168,19.2,0,24,0) ; z = -z0 =1730.4; lag.ineqlin =(1.58;3.26; 0.00) ; … shili0801.m
实例1:奶制品生产销售计划 实例1:奶制品生产销售计划 Max==12x1+8x2+22x3+16x4-1.5x5-1.5 Max二=12x+8x2+22x2+16x-1.5x-1.5x ≤50 曰4x1+3x2 +4xs+3x6≤ 五+五+≤50 4x1+x3)+2(x2+x6)+2x+2x6≤480 4(x1+x)+2(x2+x)+2x,+2x,≤480口2x+x2 +3xs+2x6≤24 x1+x5≤100 zl=1731.98 1-z=1731.98-17304=1.58 x,x2,其3,x4,x5,x620 x1,x2,x3,x4,x5,x620 所以1小时劳动时间的影 子价格应为3.26/2=1.6 0,168,19.2,0,24,0);==-20=1730.4 “影子价格” =(0,168,192,0,24,0);二=-20=1 单位劳动时间增加的利 inulin=(1.58:3.26;0.00) lag. ineqlin =(1.58; 3. 26, 0.00 润是1.63(元 l*12=1.58·12=18.96>15 →15元可增加桶牛奶,应否投资? 应该投资1 聘用临时工人增加劳动时间,工资最多每小时几元? 大学学实编) 实例:奶制品生产销售计划 NLP基本原理(不等式约束) Max==12x1+8x2+22x+16x4-1.5x-1.5x6 min ==f(x) 设x为可行解 若每公斤B的获利下降10% 位于约束边界/g0 4x+x)+2x2+x)+2+2x5480应将目标函敷中x的系数改 st.g(x)≤0,j= 为198,重新计算发现最优 x+x≤100 解和最优值均发生了变化 (x)=0.j∈J起作用约束(j=1) 若B2的获利向上波动10%, g(x)<0.j∈J不起作用约束j=23 4=075 原计划也不再是最优的 可行方向dx+M∈G(0<A<) MATLAB没有给出这种敏 81(x+d)=g,(x)+Vg(x)d+O(2) x=0.1681920240:=0=1704感性分析的结果(LNDC g1(x)2d<0(∈h1)( lag. ineqlin=(1. 58: 3. 26: 0.00) 和 LINGO软件可以给出) 下降方向df(x+)<f(x)0<<4)vf(xyd<0(2) ·B1,B2的获利经常有10%的波动,对计划有无影响? 着沿方向既可行又下降,则不是最优解③aa (学学奖 大学费学买验 Te Vgj(x)d<0 (eJn()v/(r)'d<0(2 KKT条件的几何解释 x为最优解口不存在满足(1)、2)的d vf(x)+∑AVg,(x)=0 minf(x)=(x-7)+(x2-3) 最优解的必要条件若x为最优解 且Vg,(x)O∈J)线性无关,则存在A1,…元1≥0 2,g,(x)=0,j=12 stg(x)=x2+x2-10≤0 8(x)=x+x2-4≤0 vf(x)+∑2Vg,(x)=0 KT条件 1g,(x)=0,j=12… 互补性条件 (7,3)。最优解在P(3,1)取得 v,Vg/(∈小)线性无关, s.b1(x)=0 Vf(P)+Vg, (P)+2Vg,(P)=0 8(x)≤0,/=1-1则存在和方20 P1)是KKT点 V(x)+Vb(x)+gvg;(x)=04(x)=0j=1…1 x其它点如Q均不是
5 • 15元可增加1桶牛奶，应否投资？ 实例1: 奶制品生产销售计划 1 2 3 4 5 6 Max z = 12x + 8x + 22x +16x −1.5x −1.5x 50 3 4 1 5 2 6 ≤ + + x + x x x 4( ) 2( ) 2 2 480 x1 + x5 + x2 + x6 + x5 + x6 ≤ x1 + x5 ≤ 100 3 5 x = 0.8x 4 6 x = 0.75x 0 x1 ,x2 ,x3 ,x4 ,x5 ,x6 ≥ x=(0,168,19.2,0,24,0) ; z = -z0 =1730.4 lag.ineqlin =(1.58;3.26; 0.00) ; … 4x1 + 3x2 + 4x5 + 3x6 ≤ 600 601 z1=1731.98 z1-z = 1731.98-1730.4=1.58 z1=lag.ineqlin(1) z1*12=1.58*12= 18.96>15 应该投资！ “影子价格” 实例1: 奶制品生产销售计划 • 聘用临时工人增加劳动时间，工资最多每小时几元？ 1 2 3 4 5 6 Max z = 12x + 8x + 22x +16x −1.5x −1.5x 50 3 4 1 5 2 6 ≤ + + x + x x x 4( ) 2( ) 2 2 480 x1 + x5 + x2 + x6 + x5 + x6 ≤ x1 + x5 ≤ 100 3 5 x = 0.8x 4 6 x = 0.75x 0 x1 ,x2 ,x3 ,x4 ,x5 ,x6 ≥ x=(0,168,19.2,0,24,0) ; z = -z0 =1730.4 lag.ineqlin =(1.58;3.26; 0.00) ; … 4x1 + 3x2 + 4x5 + 3x6 ≤ 600 2 3 2 240 x1 + x2 + x5 + x6 ≤ lag.ineqlin(2)=3.26， 所以1小时劳动时间的影 子价格应为3.26/2=1.63， 即单位劳动时间增加的利 润是1.63(元) • B1，B2的获利经常有10%的波动，对计划有无影响？ 实例1: 奶制品生产销售计划 1 2 3 4 5 6 Max z = 12x + 8x + 22x +16x −1.5x −1.5x 50 3 4 1 5 2 6 ≤ + + x + x x x 4(x1 + x5 ) + 2(x2 + x6 ) + 2x5 + 2x6 ≤ 480 x1 + x5 ≤ 100 3 5 x = 0.8x 4 0 75 6 x = . x 0 x1 ,x2 ,x3 ,x4 ,x5 ,x6 ≥ 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%，y 原计划也不再是最优的 MATLAB没有给出这种敏 感性分析的结果（LINDO 和LINGO软件可以给出） NLP基本原理(不等式约束) st g x j l z f x j x n . . ( ) 0, 1,..., min ( ) ≤ = = ∈ℜ ∇g1 d g2=0 g1=0 g3=0 gj <0 o x 设x为可行解， 位于约束边界 1 g (x) 0, j J j = ∈ J1~起作用约束(j=1) 2 g (x) 0, j J j < ∈ J2~不起作用约束(j=2,3) (0 ) +λ ∈G <λ <λ0 可行方向d x d 下降方向d ( ) ( ) (0 ) +λ < < λ < λ0 f x d f x (G) ( ) 0 ( ) (1) 1 g x d j J T ∇ j < ∈ ∇f (x) d < 0 (2) T ( ) ( ) ( ) ( ) 2 g x λd g x λ g x d O λ T j + = j + ∇ j + 若x沿d方向既可行又下降，则x不是最优解 最优解的必要条件 ( ) ( ) 0 1 ∇ + ∑ ∇ = = f x g x j l j λi g x j l λ j j( ) = 0, =1,2,L ~KKT条件 互补性条件 g x j l s t h x i m f x j i x n ( ) 0 , 1,..., . . ( ) 0 , 1,..., min ( ) ≤ = = = ∈ℜ 则存在 和 ， 线性无关， 0 , ( ) 1 ≥ ∇ ∇ ∈ i j i j h g j J µ λ ( ) ( ) ( ) 0 1 1 ∇ + ∑ ∇ + ∑ ∇ = = = f x h x g x j l j i m i µ i i λ g x j l λ j j( ) = 0, =1,L x为最优解 不存在满足(1),(2)的d ( ) 0 ( ) (1) 1 g x d j J T ∇ j < ∈ ∇f (x) d < 0 (2) T 且 线性无关，则存在 ∇g j (x) ( j ∈ J 1 ) λ1 ,Lλ l ≥ 0 若x为最优解， KKT条件的几何解释 Q P ∇f -∇g1 -∇g2 ∇f -∇g1 -∇g2 ( ) 0 ( ) 4 0 . . ( ) 10 0 min ( ) ( 7) ( 3) 3 2 2 1 2 2 2 2 1 1 2 2 2 1 = ≥ = + − ≤ = + − ≤ = − + − g x x g x x x st g x x x f x x x ( ) ( ) 0 1 ∇ + ∑ ∇ = = f x g x j l j λi g x j l λ j j( ) = 0, =1,2,L x2 x 0 1 最优解在P(3,1)取得 ( ) ( ) 2 ( ) 0 1, 2, 0, 1 2 1 2 3 ∇ +∇ + ∇ = = = = f P g P g P λ λ λ P(3,1)是KKT点 其它点(如Q)均不是 • (7,3)
二次规划(QP)及有效集方法 解二次规划的有效集方法 minf(x)=xHx+cx当H为对称阵,称二次规划 基本思想:对于不等式约束的二次规划,在某可行点 当H正定时,称凸二次规划 处将不起作用约束去掉,起作用约束视为等式约束 st.Ax≤b 凸二次规划性质: 通过求解等式约束的二次规划来改进可行点 等式约束Ax=b下 最优点←KKT点 的 Lagrange乘子法部最优解全局最优解 min(x)=-r'Hx+cx minf(x)=5x'Hx+ s.Ax≤b L函数L(x,A)=xHx+ax+(Ax-b),A∈R b,j∈J(a,是A的第/列) 基本·若x为(1)的最优解,则它也是(2)的最优解 最优解方程x+c+4=0 原理若x为(1)的可行解,又是(2)的KT点,且 Ax-b=0 L一乘子非负,则它必是(1)的最优解 基本步骤设(1)的可行点为x*,有效集 min==lxX+cx 记作/*,用L—乘子法求解 MTAB求解QPs,Ax≤b,.x=b:5xn Ix, fval, exitflag, output, lambda] minf(x*d)=1(x*+d)H(x*+d)+(x*+4)得P”, quadprog(H,C, Al, bl, A2, b2, vl, v2, xO, options) 若d=0,则x*为(2)最优解,当λ*非负时x*是(1)最优解 x1+2x2=3 x1-x22-3H 有效集若=0,且(2)<0,q∈/,则x不是最优解, x1-3x,≤3 修正有效集修正为g} x1≥2,x2≤04 若d≠0,以此为方向确定步长a·使得 a(x*+*弹)=by,pgJ*,则有效集修正为严U{p Exam0803 m (学学奖 实例2:投资组合问题 实例2:投资组合问题 minz2=4x2+36x2+100x2+5x1x2-20xx3-30x2x1 加权 MinZ=BZ2-Z st5x1+8x2+10x321000 0x1+25x2+30x3≤5000 模型 st.20x1+25x2+30x3≤5000 portfolio2.m 通过试探发现B从 解得x=1.0e+002*(13111,0.1529,0.2221) 0.0001-01以0.0001的 如果一定要整数解,可以四舍五入到(131,15,22) 步长变化就可以得到 如利用 LINGO软件可得整数最优解(132,15,2) 很好的近似结果 用去资金为132x20+15×25+22×30=3675(百元) 期望收益为132×5+15×8+22×10=1000(百元 投资股纂A、B、C分别 风险(方差为68116,标准差约为261(百元 为153、35、35(手)
6 二次规划(QP)及有效集方法 st Ax b f x x Hx cx T ≤ = + . . 2 1 min ( ) 当H为对称阵，称二次规划 当H正定时，称凸二次规划 凸二次规划性质： 最优点⇔KKT点； 局部最优解⇔全局最优解； T T m L x λ = x Hx + cx + λ (Ax − b), λ ∈ R 2 1 ( , ) 最优解方程 0 0 − = + + = Ax b Hx c A T T λ L函数 等式约束 下 的Lagrange乘子法 Ax = b 解二次规划的有效集方法 基本思想：对于不等式约束的二次规划，在某可行点 处将不起作用约束去掉，起作用约束视为等式约束， 通过求解等式约束的二次规划来改进可行点。 • 若x为(1)的最优解，则它也是(2)的最优解 • 若x为(1)的可行解，又是(2)的KKT点，且 L—乘子非负，则它必是(1)的最优解。 st Ax b f x x Hx cx T ≤ = + . . 2 1 min ( ) (1) 1 . . , 2 1 min ( ) st a x b j J f x x Hx cx j j T = ∈ = + (a 是A的第 j列) j (2) 基本 原理 . . 0, * ( * ) ( * ) ( * ) 2 1 min ( * ) st a d j J f x d x d H x d c x d j T = ∈ + = + + + + 若d*≠0，以此为方向确定步长α*使得 ap(x*+α*d*)=bp，p∉J*，则有效集修正为J*∪{p}。 设(1)的可行点为x*，有效集 记作J*，用L—乘子法求解: 基本步骤 得d*, λ* • 若d*=0, 则x*为(2)最优解; 当λ *非负时x*是(1)最优解 若d*=0, 且(λ *)q<0, q∈J* , 则x*不是最优解, 有效集修正为J*\{q} 。 有效集 修正 MATLAB求解QP [x,fval,exitflag,output,lambda] = quadprog(H,c,A1,b1,A2,b2,v1,v2,x0,options) 1 1 2 2 1 2 2 1 . . , , min s t A x b A x b v x v z x Hx c x T T ≤ = ≤ ≤ = + 2, 0 3 3 2 3 . . 2 3 min ( , ) 2 3 3 3 1 2 1 2 1 2 1 2 1 2 2 1 2 2 2 1 2 1 ≥ ≤ − ≤ − ≥ − + = = − + − + x x x x x x st x x f x x x x x x x x Exam0803.m (1,2), 3, [2, Inf], 2 [Inf,0] , 3 3 , 1 3 2 1 1 3 , 3 6 4 3 2 2 1 1 1 = = = − = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛− = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = A b v v A b H c 解得x = 1.0e+002 *（1.3111，0.1529，0.2221） 如果一定要整数解，可以四舍五入到（131，15，22） 如利用LINGO软件,可得整数最优解(132，15，22) 用去资金为132×20+15×25+22×30 = 3675（百元） 期望收益为132×5+15×8+22×10 = 1000（百元） 风险(方差)为 68116，标准差约为261（百元） 实例2：投资组合问题 s.t. 5x1 +8x2+10x3 ≥ 1000 20x1+25x2+30x3 ≤ 5000 x1，x2，x3 ≥ 0 1 2 1 3 2 3 2 3 2 2 2 2 1 min Z = 4x + 36x +100x + 5x x − 20x x − 30x x portfolio1.m 通过试探发现 β从 0.0001~0.1以0.0001的 步长变化就可以得到 很好的近似结果 实例2：投资组合问题 Min Z =βZ2 - Z1 s.t. 20x1+25x2+30x3 ≤ 5000 x1，x2，x3 ≥ 0 加权 模型 0 200 400 600 800 1000 1200 1400 1600 1800 0 100 200 300 400 500 600 700 800 900 预期收益(百元) 均方差（百元） 投资股票A、B、C分别 为153、35、35（手） portfolio2.m
非线性规划LP)的解法 可行方向法、罚函数法、梯度投影法 min d 'Gd+v(x,)'d 逐步二次规划法( Sequential Quadratic Programming) 题 s1. Vh(x,d+h(x n f(x) vg,(x,)d+g,(x)≤0,j=1…1 SQP的基本原理 x是第k次选代的初始点,G是海赛阵VL的近似 构造LNP的拉格朗日函数 gj(x)≤0,j=1 将最优解d作为选代的搜索方向,令x+1=X+akd L(x,,A)=f(x)+2h(x)+∑g,(x) SOP的·求解QP子问题,得d; 用二次函数近似L(x,A)NLP化为QP 基本步骤·用线性搜索计算迭代步长a 再解一系列QP子问题 ·确定矩阵G的迭代公式 min ==f(r) MATLAB求解 =f(x) MATLAB求解s.c1(x)≤0,c2(x)=0, 约束MLP s1.c1(x)≤0,c2(x)=0, 约束MLP A1x≤b1,A2x=b2,V1≤x≤V2 A1x≤h,A2x=b2,n≤x≤v2 alcon.m给出约束, Gradconstr='on时还给出梯度,形式为 Ix, fval, exitflag, output, lambda, grad, hessian] function [cl, c2, GCl, GC2]= alcon(x) rincon(@fun, xO, Al, bl, A 2, b2, vl, v2, @nIcon, options, Pl, P2,. 8 nonlinear inequalities at x fun.m给出函数f,当 Gradobj='on’时必须给出f的梯度 8 nonlinear equalities at x Hessian='on’时还必须给出其 Jacobi矩阵,形式为 if nargout 2 =.. 8 gradients of cl function [f, g, H] fun(x) s gradients of c2 8 objective function value if nargout >1 2x2+1) x2-x2+1.5≤0,xx2+1020 x2+x2-1=0 Exam0s04, m (学学奖 大学费学买验 实例3:供应与选址三 实例3:供应与选址 决策变量tc 分之x-a)2+(-4)32 运量)-12维4∑-4,1=1.6 S之9xy-a2+-)}2改建两个新料场 决策变量 cy=d,i=1-6 c,(xy)-16维 线性规划模型 Shili0803. m: 用例中数 cy se, j=1,2 据计算 c,(料场A) 初始点的选择 最优解为c。(料场B) 用现料场总吨公里数为1362 局部最优解问题 总吨公里数为1362 Shilioso3linn 结果:总吨公里数为853,比使用原料场减少了50.9
7 非线性规划(NLP)的解法 可行方向法、罚函数法、梯度投影法... 逐步二次规划法(Sequential Quadratic Programming) SQP的基本原理 构造LNP的拉格朗日函数 ( , , ) ( ) ( ) ( ) 1 1 L x f x h x g x j l j j m i ∑ i i ∑ = = µ λ = + µ + λ g x j l s t h x i m f x j i x n ( ) 0, 1,..., . . ( ) 0, 1,..., min ( ) ≤ = = = ∈ℜ 用二次函数近似 L(x,µ,λ), NLP化为QP, 再解一系列QP子问题。 QP子 问题 g x d g x j l st h x d h x i m d G d f x d j k T j k i k T i k T k k T L L ( ) ( ) 0, 1, . . ( ) ( ) 0, 1, ( ) 2 1 min ∇ + ≤ = ∇ + = = + ∇ k x 是第k 次迭代的初始点，Gk 是海赛阵 L 2 ∇ 的近似。 • 求解QP子问题，得 dk； 将最优解dk作为迭代的搜索方向，令 k k k k x = x +α d +1 SQP的 基本步骤 • 确定矩阵Gk的迭代公式。 • 用线性搜索计算迭代步长αk； MATLAB求解 约束NLP [x,fval,exitflag,output,lambda,grad,hessian] = fmincon(@fun,x0,A1,b1,A2,b2,v1,v2,@nlcon,options,P1,P2, ...) fun.m给出函数f，当GradObj=‘on’时必须给出f的梯度， Hessian=‘on’时还必须给出其Jacobi矩阵，形式为 1 1 2 2 1 2 1 2 , , . . ( ) 0, ( ) 0, min ( ) A x b A x b v x v s t c x c x z f x ≤ = ≤ ≤ ≤ = = function [f,g,H] = fun(x) f = ... % objective function value if nargout > 1 g = ... % gradient of the function if nargout > 2 H = ... % Hessian end nlcon.m给出约束,GradConstr=‘on’时还给出梯度,形式为 1 1 2 2 1 2 1 2 , , . . ( ) 0, ( ) 0, min ( ) A x b A x b v x v s t c x c x z f x ≤ = ≤ ≤ ≤ = = 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 end MATLAB求解 约束NLP 例4 Exam0804.m min (4 2 4 2 1) 1 2 2 2 2 2 1 1 e x + x + x x + x + x 1.5 0, 10 0 x1 x2 − x1 − x2 + ≤ x1 x2 + ≥ 1 0 2 2 x1 + x − = 用例中数 据计算， 最优解为 i 123456 i1 c （料场 A） 350701 i 2 c （料场 B） 0 0 4 0 6 10 总吨公里数为136.2 , 1,2 . . , 1,...,6 min [( ) ( ) ] 6 1 2 1 2 1 6 1 2 2 1/ 2 ≤ = = = − + − ∑ ∑ ∑∑ = = = = c e j s t c d i c x a y b ij j i ij i j j i ij j i j i 线性规划模型 决策变量：ci j (料场j到工地i的 运量）~12维 Shili0803lin.m 实例3: 供应与选址 结果：总吨公里数为85.3，比使用原料场减少了50.9。 初始点的选择 实例3: 供应与选址 , 1,2 . . , 1,...,6 min [( ) ( ) ] 6 1 2 1 2 1 6 1 2 2 1/ 2 ≤ = = = − + − ∑ ∑ ∑∑ = = = = c e j st c d i c x a y b ij j i ij i j j i ij j i j i 决策变量： ci j，(xj ,yj )~16维 用现料场总吨公里数为136.2 改建两个新料场 局部最优解问题 Shili0803.m; shili083fun.m
学学实纷 计算方法的改善 新料场位置〔x,y) min (c(n-ai)2+(1-b)21/2 05000 (7.24977.7499) (d1-c1(x2-a)2+(y2-b)2y2 ≤ ≤q,)(d1-c)≤e2 决簟变量:c,(xy}-10维 局部最优解问题有所改进 +为工地,数字为用量,*为新料场,数字为供应量。 布置实验内容 1)掌握用 MATLAB优化工具包解线性规划和 实验非线性规划(包括二次规划的方法 目的」2)练习建立实际问题的线性规划和非线性 规划模型 内容课上布置,或参见网络学堂 还 8
8 2 6 1 1 6 1 2 1/ 2 2 2 2 6 1 2 1/ 2 1 2 1 , ( ) . . , 1, 6 ( )[( ) ( ) ] } min { [( ) ( ) ] c e d c e s t c d i d c x a y b c x a y b i i i i i i i i i i i i i i i ≤ − ≤ ≤ = + − − + − − + − ∑ ∑ ∑ = = = L 决策变量：ci ，(xj ,yj )~10维 计算方法的改善 局部最优解问题有所改进 Shili0803n.m; shili083fun1.m i 12345 6 新料场位置（ , ) j j x y i1 c 30476 0 （3.2552 5.6528） i 2 c 0 5 0 0 0 11 （7.2497 7.7499） +为工地, 数字为用量; *为新料场, 数字为供应量。 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 3 5 4 7 6 11 20 16 实验 目的 1）掌握用MATLAB优化工具包解线性规划和 非线性规划(包括二次规划)的方法； 2）练习建立实际问题的线性规划和非线性 规划模型。 实验 内容 课上布置，或参见网络学堂 布置实验内容