《经济数学基础》课程教学资源：实验7 无约束优化.pdf_大学文库

2 1 2 1 1 1 2 2 2 , max ( , ) ( ) ( ) 1 2 z x x p q x p q x x x = − + − 目标利润最大 , , , 0 , , , 0 2 2 2 2 2 2 1 1 1 1 1 1 2 2 1 1 = + > = + > − − q r e c r c q r e c r c x x λ λ λ λ 假设C q随x (本牌号)增加而减小，呈负指数关系实例1 产销量安排 0 y x VOR2 x=629, y=375 309.00 (1.30) 864.3(2.0) 飞机 x=?, y=? VOR1 x=764, y=1393 161.20 (0.80) VOR3 x=1571, y=259 45.10 (0.60) 北 DME x=155, y=987 飞机与监控台（图中坐标和测量距离的单位是“公里”）实例2 飞机精确定位问题飞机精确定位模型 4 2 4 2 4 ( ) ( ) atan2( , ) x x y y d x x y y i i i − + − = − − = θ 不考虑误差因素超定方程组，非线性最小二乘！要求计算：飞机位置坐标（）记测量距离为，距离误差为记测量角度为，角度误差为已知数据：设备位置坐标分别为 x y d i x y i i i i i , . , 1,...,3; ( , ), 1,...4; 4 σ 4 θ σ = = 量纲不符！？ 2 4 2 4 2 4 3 1 2 atan2( , ) ( ) ( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + − + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − ∑= d x x y y x x y y Min i i i i x,y θ 飞机精确定位模型 4 4 2 4 2 4 4 4 ( ) ( ) atan2( , ) σ σ θ σ θ σ − ≤ − + − ≤ + − ≤ − − ≤ + d x x y y d x x y y i i i i i i 考虑误差因素 Min x; Min y; Max x; Max y. 以距离为约束，优化角度误差之和（或平方和）；或以角度为约束，优化距离误差. 非线性规划误差非均匀分布！？？仅部分考虑误差! 角度与距离的“地位”不应不同！有人也可能会采用其他目标，如：飞机精确定位模型 2 4 2 4 2 4 4 2 3 1 atan2( , ) ( ) ( ) ( , ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − + − + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − − = ∑= σ σ x x y y θ d x x y y Min E x y i i i i i 误差一般服从什么分布？正态分布！不同的量纲如何处理？无约束非线性最小二乘模型归一化处理！角度需要进行预处理，如利用atan2函数, 值域(-pi, pi) 无约束优化:最优解的分类和条件 Min f (x) x 给定一个函数 f(x),寻找 x* 使得 f(x*)最小，即 T n n 其中 x = ( x1 , x 2 ,L , x ) ∈ ℜ 局部最优解全局最优解必要条件 ( ) ( , , ) 0 1 * ∇ = = T x xn f x f L f x * f(x) xl xgo 充分条件 ( ) 0, ( ) 0 * 2 * ∇f x = ∇ f x > Hessian阵 n n i j x x f f × ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∇ = 2 2

3 求解无约束优化的基本思路基本思想在 n ℜ 中某一点，确定一个搜索方向及沿该方向的移动步长，得到使目标函数下降的新的点迭代步骤 Step 1 初始化：初始点x0，终止准则等 Step 2 迭代改进：方向d k，步长α k ( ) ( ) k 1 k f x < f x k k k k + x = x + α d +1 Step 3 终止检验：得到近优解或k+1⇒k转2 选择d k ,α k 使 f 下降更快 ⇒ 不同算法搜索方向的选择 1 最速下降法（梯度法）下降方向最速下降方向迭代改进格式算法特点初始阶段改进较快，最优解附近改进较慢 k k k k T k k f (x ) f (x d ) f (x ) f (x )d 1 = + = + ∇ + ∇ ( ) < 0 T k k f x d ( ) k k d = −∇f x ( ) k 1 k k x = x − ∇f x + 将 ( ) k +1 f x 在 k x 点作泰勒展开，只保留一阶项，有暂不考虑搜索步长，可设αk=1 (负梯度方向) 2 Newton方法特点局部2阶收敛;需计算Hessian阵,它可能病态或不正定将f(xk+1)在xk点作泰勒展开至二阶项，用d替代dk 牛顿方程牛顿方向迭代格式 ( ) ( ) 2 k k ∇ f x d = −∇f x f x f x d f x f x d d f x d k k k T k T k ( ) 2 1 ( ) ( ) ( ) ( ) 1 2 = + = + ∇ + ∇ + ( ( )) ( ) k 2 k 1 k d = − ∇ f x ∇f x − ( ( )) ( ) k 1 k 2 k 1 k x = x − ∇ f x ∇f x + − [ ( )] ( ) k 1 k k 1 k x x F x F x + − 比较 = − ′ 的牛顿法解F ( x) = 0 F ( x) ⇔ ∇f ( x) 求d使f(xk+1)极小⇒右端对d导数为0 ( ) ( ) 0 2 ⇒∇f x +∇ f x d = k k 3 拟Newton方法目的不计算Hessian阵，克服病态、不正定、计算复杂等缺陷，同时保持收敛较快的优点回顾解方程组 F(x)=0的拟牛顿法 ( ) ( ) ( ) −1 −1 − = − k k k k k k 使A 满足 A x x F x F x k k k k 用迭代方法A A A 计算A −1 −1 = + ∆ 思路 [ ( )] ( ) k 1 k k 1 k x x F x F x + − = − ′ ( ) ( ) k 1 k k 1 k x x A F x + − = − Min f (x) 优化问题 x ∇f ( x)相当 F ( x) f x F x f G f 2 2 2 ∇ ( )相当 ′( ), ∇ 不一定正定，构造正定阵代替∇ 3 拟Newton方法(续) k k k k k k G = G + ∆G H = H + ∆H 构造迭代公式 +1 或 +1 按照拟牛顿条件： k k k k k k G ∆x = ∆f ∆x = H ∆f +1 或 +1 , ( ) ( ) k k 1 k k k 1 k ∆x = x − x ∆f = ∇f x − ∇f x 记 + + ( ) k 1 k k k x = x − H ∇f x + 设在第k步, Gk已得到, Hk=(Gk)-1,可计算 ( ) +2 +1 +1 +1 = − ∇ k k k k 于是有 x x H f x 3.1 Davidon-Fletcher-Powell(DFP)公式 k T k k k k k T k k T k k k T k f H f H f f H x f x x H ∆ ∆ ∆ ∆ − ∆ ∆ ∆ ∆ ∆ = ( ) ( ) ( ) ( ) 3 拟Newton方法(续) k T k k k T k k k k T k T k k k T k T k k T k k k x f f x G G x f x f f f x f x G x G ∆ ∆ ∆ ∆ + ∆ ∆ − ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ = + ( ) ( ) ( ) ( ) ( ) ) ( ) ( ) (1 3.2 Broyden-Fletcher-Goldfarb-Shanno(BFGS)公式 k T k k k k k T k k T k k k T k x G x G x x G f x f f G ∆ ∆ ∆ ∆ − ∆ ∆ ∆ ∆ ∆ = ( ) ( ) ( ) ( ) k T k k k T k k k k T k T k k k T k T k k T k k k x f x f H H f x x f x x x f f H f H ∆ ∆ ∆ ∆ + ∆ ∆ − ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ = + ( ) ( ) ( ) ( ) ( ) ) ( ) ( ) (1

4 搜索步长的确定线性搜索(Line Search)确定步长方法问题给定xk和方向dk, 确定步长αk, 使得 min ( ) k k f x αd α + 优化算法黄金分割(0.618)法、Fibonacci法、Newton 切线法、割线法、2次或3次插值法等 ~一维优化问题 k = 0 d df α α ( ) ∇ ( + ) = 0 k T k k k d f x α d 非线性最小二乘(Least Square)拟合问题给定(t i , y i ), i=1,…n, 拟合一个函数y=f(t, x)，其中x为待定的参数向量, f 对x非线性。优化模型 min R(x) r (x)r(x) T x = r(x) y f (t , x) i = i − i T n r(x) (r(x), r (x)) 记误差 = 1 L 根据目标函数是 r(x) 的二次函数的特点构造简单算法非线性最小二乘拟合 R 2J (x) r(x) T ∇ = R J x J x S T 2 ( ) ( ) 2 2 ∇ = + i i k l m m r x r x x × ∇ ( ) = (∂ ∂ ∂ ) 2 2 讨论 • 牛顿法要计算Hessian矩阵，其中S计算量大 • 若 f 对x线性, 则化为线性最小二乘拟合,此时S=0 特定算法考虑如何忽略或近似矩阵S。 i j n m J x r x × 记 r(x) 的雅各比阵为 ( ) = (∂ / ∂ ) R(x) r (x)r(x) T = ( ) ( ) 2 1 S r x r x i n i =∑ i ∇ = Gauss-Newton算法：忽略矩阵S ( ) ( ) ( ) ( ) k T k k k k J x J x d = −J x r x f用R代替，下降方向dk满足 G-N算法收敛性依赖 f 对 x 的线性程度, 及偏差r的大小非线性最小二乘拟合 ( ) ( ) 2 k k 牛顿方程 ∇ f x d = −∇f x R 2J (x) r(x) T ∇ = 2 ( ) ( ) 2R J x J x T ∇ = Levenbery-Marquardt算法：G-N算法修正 ( ( ) ( ) ) ( ) ( ) k T k k k k k J x J x +α I d = −J x r x 其中αk>0为修正参数. dk位于牛顿方向(αk很小)和负梯度方向(αk很大)之间 L-M算法非线性最小二乘拟合 ( ) ( ) ( ) ( ) k T k k k k J x J x d = −J x r x 防止 JTJ 出现病态优化工具箱基本用法： x=fminbnd(@f,lb,ub) x=fminunc(@f,x0) x=fminunc(@f,x0,options,P1,P2,...) x=fminsearch(@f,x0,options,P1,P2,...) f.m ~ f(x)的m文件名 x0 ~初始点; x ~最优解 P1,P2, … ~ 传给fun的参数中间输入项缺省用[ ]占位 n x 模型：Min f (x), x∈ R 2 2 min 2 2 = = + a b b y a x 其中例： Exam0701.m Exam0702.m x [1,8] 1 min (3sinx x) ∈ 例： +

5 非线性最小二乘法 min R(x) r (x)r(x) T lb x ub = ≤ ≤ r ( x) y f (t , x) i = i − i T n r(x) (r (x), r (x)) = 1 L [x,resnorm,res,exitf,out,lambda,jacob]= lsqnonlin(@fun,x0,lb,ub,options,P1,P2,…) 输入的用法与fminunc类似，但注意： f.m~f(x)的m文件名: function y=f(x,t) fun.m~r(x)的m文件名: function r=fun(x,t,y) 输出 resnom=r(x)T*r(x), res=r(x)(误差向量) [x,resnorm,res,exitf,out,lambda,jacob]= lsqcurvefit(@f,x0,t,y,lb,ub,options,P1,P2,…) 非线性最小二乘法 min R(x) r (x)r(x) T lb x ub = ≤ ≤ c 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01 t 0.25 0.5 1 1.5 2 3 4 6 8 kt c t re − ( ) = 例3 用下面一组数据拟合系数r, k : exam0703fit.m exam0703lsq.m 控制精度，观察中间结果，控制迭代次数等 Optimset //显示控制参数 optimset optfun //显示optfun的控制参数 opt=optimset //控制参数设为[](即缺省值) opt=optimset(optfun)//optfun控制参数设缺省值 Opt=optimset('par1',val1,'par2',val2,...) Opt=optimset(oldopts,'par1',val1,...) opt=optimset(oldopts,newopts) val=optimget(opt,'par1','par2',…) val=optimget(opt,'par1','par2',…, default) Options 控制参数设定/获取: optimset; optimget Diagnostics ‘on’ | {‘off’} //是否显示诊断信息 Display 'off' | 'iter' | ‘final’ | ‘notify‘ //显示信息的级别 GradObj ‘on’ | {‘off’}//是否采用分析梯度 Jacobian ‘on’ | {‘off’} //采用分析Jacob阵（用于约束优化中） LargeScale ‘on’ | {‘off’}//是否采用大规模算法 MaxFunEvals 最大函数调用次数 MaxIter 最大迭代次数 TolCon 约束的控制精度（用于约束优化中） TolFun 函数值的控制精度 TolX 解的控制精度主要控制参数(对大规模/中等规模算法均有效) 最一般的输出形式 [x,f,exitflag,out,grad,hess]=fminunc(...) 更多输出：最优值等 f 目标函数值 exitflag >0收敛,0达到函数或迭代次数, <0不收敛 Output iterations 实际迭代次数 funcCount 实际函数调用次数 algorithm 实际采用的算法 cgiterations 实际PCG迭代次数（大规模算法用） stepsize 最后迭代步长（中等规模算法用） firstorderopt 一阶最优条件（梯度的范数） grad 目标函数的梯度 hess 目标函数的Hessian矩阵例 4 min , 10, 1 2 2 2 1 + a = b = b x a x Exam0704.m 用optimset控制参数选择

6 算法选择：optimset中参数控制 LargeScale ‘on’ | ‘off’ （on为缺省）搜索步长的算法选择系统缺省采用大规模算法（如果可能）当设为‘off’时缺省：BFGS、混合2, 3次多项式插值 LineSearchType='quadcubic' 混合2, 3次多项式插值 LineSearchType='cubicpoly' 3次多项式插值 HessUpdate = ‘dfp’ （DFP算法） HessUpdate = ‘steepdesc’（最速下降算法）搜索方向的算法选择计算结果方向步长最优解x 最优值f n BFGS 2,3次 (9.9997e-001 9.9994e-001) 1.0944e-009 165 DFP 2,3次 (9.9997e-001 9.9994e-001) 2.1173e-009 165 steep 2,3次 (-8.0263e-001 6.5064e-001) 3.2536e+000 1003 BFGS 3次 (1.0001e+000 1.0003e+000) 3.1432e-008 162 DFP 3次 (9.0468e-001 8.1572e-001) 9.8270e-003 495 steep 3次 (-1.1831e+000 1.4056e+000) 4.7695e+000 1002 2 2 2 例5. min f (x, y) =100(y − x ) + (1− x) 精确解：x=y=1, f(x,y)=0 Exam0705.m 无约束优化 n x 模型：Min f (x), x∈ R 采用分析梯度： GradObj='on' x=fminunc(@fun,x0,opt,…) fun.m中还要有一般形式 function [f,g]=fun(x) ∇f ( x) 2 2 2 续例5. min f (x, y) =100( y − x ) + (1− x) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − − − − − ∇ = 200( ) 400 ( ) 2(1 ) 2 2 y x x y x x 算出 f exam0705grad_run.m 方向步长最优解x 最优值f n BFGS 2,3次 (1.0001e+000 1.0002e+000) 8.9857-009 105 DFP 2,3次 (1.0001e+000 1.0003e+000) 2.2499e-008 109 BFGS 3次 (1.0000e+000 1.0001e+000) 2.7164e-009 53 DFP 3次 (7.4711e-001 5.5212e-001) 6.7609e-002 144 计算结果与不用分析梯度的结果比较方向步长最优解x 最优值f n BFGS 2,3次 (9.9997e-001 9.9994e-001) 1.0944e-009 165 DFP 2,3次 (9.9997e-001 9.9994e-001) 2.1173e-009 165 BFGS 3次 (1.0001e+000 1.0003e+000) 3.1432e-008 162 DFP 3次 (9.0468e-001 8.1572e-001) 9.8270e-003 495 x0 * O f (x, y) =100(y − x2 ) 2 + (1− x) 2 的图形 exam0705plot.m Rosenbrock 函数（香蕉函数）算法选择：BFGS公式，混合2，3次插值，一般较好。几个值得注意的问题精度控制：对迭代次数有重大影响，应适当选择。梯度函数：利用分析梯度可能改进算法的性能改变初始值由一个初值出发通常得到局部最优解，如果函数存在多个局部最优，只有改变初值，对局部最优进行比较，才有可能得到全局最优解。其他算法选择：（详细用法请查阅help文档）高度非线性、不连续时可用程序 fminsearch(@fun,x0) 单变量时可用程序 fminbnd(@fun,v1,v2)

7 非线性最小二乘法 min R(x) r (x)r(x) T lb x ub = ≤ ≤ r ( x) y f (t , x) i = i − i T n r(x) (r (x), r (x)) = 1 L [x,resnorm,res,exitf,out,lambda,jacob]= lsqnonlin(@fun,x0,lb,ub,options,P1,P2,…) 输入的用法与fminunc类似，但注意：fun.m ~r(x) 的m文件名，Jacobian=‘on’时含有导数信息 function [r,J] = fun(x) ∇r(x) 算法选择: 缺省：大规模算法（LargeScale = 'on' ）当LargeScale = 'off' ： •缺省： Levenberg-Marquardt算法； •LevenbergMarquardt=‘off’:Gauss-Newton法一维搜索（线搜索）步长选择与fminunc中类似非线性最小二乘法 min R(x) r (x)r(x) T lb x ub = ≤ ≤ exam0705lsq.m exam0705lsq_jacob_run.m 实例1 产销量安排 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ =⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ 280 100 0.2 2 1 0.1 21 22 11 12 b a a a a 已知数据 = (30 100 ), = (0.015 0.02 ), = (20 30) T T T r λ c 2 21 1 22 2 2 2 2 1 11 1 12 2 1 1 1 ( ) min ( ) ( ) 2 2 1 1 b a x a x r e c x f x b a x a x r e c x x x − − − − − = − − − − − − − λ λ , , , 0, 1,2, pi =bi −ai1x1 −ai2x2 bi ai1 ai2 > i = = + , , , > 0, =1,2, − q r e c r c i i i i i x i i i i λ λ 1 2 1 1 1 2 2 2 , max ( , ) ( ) ( ) 1 2 z x x p q x p q x x x = − + − 原问题 x = 23.9025 62.4977 y = 6.4135e+003 即甲产量为23.9025，乙产量为62.4977，最大利润为6413.5 命令和最优解 x=fminunc(@f, x0) shili0701.m 其解为x1 = b1 / 2a11 = 50， x2 = b2 / 2a22 = 70 初始点的选择实例1 产销量安排 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ =⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ 0.2 2 1 0.1 21 22 11 12 a a 忽略成本及价 a a 格中的 a12, a21 1 11 1 1 2 22 2 2 问题简化为 min f (x) = −(b − a x )x − (b − a x )x 作为原问题的初始点飞机坐标(978.31，723.98), 误差平方和0.6685 (<< 4) 实例2 飞机精确定位模型 2 4 2 4 2 4 4 2 3 1 atan2( , ) ( ) ( ) ( , ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − + − + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − − = ∑= σ σ x x y y θ d x x y y Min E x y i i i i i 无约束非线性最小二乘模型角度需要进行预处理，如利用atan2函数, 值域(-pi, pi) shili0702.m 其他:非线性方程(组) 单变量方程求根： f (x) = 0, x ∈ R [x,fval,exitflag,output] = fzero(fun,x0,options,P1,P2,...) [x,fval,exitflag,output,jacobian] = fsolve(fun,x0,options,P1,P2, ... ) n 多变量方程组求零点： F(x) = 0, x ∈ R