1 大学数学实验 Experiments in Mathematics 实验7 无约束优化 清华大学数学科学系 最优化是工程技术、经济管理、科学研究、社 会生活中经常遇到的问题 , 如 : 优化模型和算法的重要意义 结构设计 资源分配 生产计划 运输方案 解决优化问题的手段 ? 经验积累，主观判断 ? 作试验，比优劣 ? 建立数学模型，求解最优策略 最优化: 在一定条件下，寻求使目标最大(小)的决策 运筹学 (OR: Operations/Operational Research) 管理科学 (MS: Management Science) 决策科学 (DS: Decision Science) (最 )优化理论是运筹学的重要内容 无 约 束 优 化 OR/ MS/ DS 优化 (Optimization), 规划 (Programming) 线 性 规 划 非 线 性 规 划 网 络 优 化 组 合 优 化 整 数 规 划 不 确 定 规 划 多 目 标 规 划 目 标 规 划 动 态 规 划 优化问题的数学模型 可行域目标函数，决策变量， ~~~ ?fx n x Rxxfz ??∈= ),(min ? 可行解（只满足 (2)）与最优解（满足 (1),(2)） ? 无约束优化（只有 (1)）与约束优化（ (1),(2)） )2(,2,1,0)(.. )1()(min mixgts xfz i x L=≤ = ? 实际问题一般总有约束，何时可用无约束优化处理？ 4. 优化工具箱的使用 2. 无约束优化的基本方法： 梯度法，牛顿法，拟牛顿法 1. 优化问题的最优解条件；算法模式 无约束优化的主要内容 3. 非线性最小二乘法 5. 实际问题中的无约束优化模型 实例1 产销量安排 牌号 产量 成本 价格 甲 x1 q1 p1 乙 x2 q2 p2 假设 A 产销平衡 假设 B p随 x (两种牌号 )增加而减小，呈线性关系 12111211121211111 ,0,,, aaaabxaxabp >>??= 某厂生产两个牌号的同一种产品，如何确定产量使利润最大 21222221222212122 ,0,,, aaaabxaxabp >>??= 2 22211121 , )()(),(max 21 xqpxqpxxz xx ?+?= 目标 利润最大 0,,, 0,,, 222222 111111 22 11 >+= >+= ? ? crcerq crcerq x x λ λ λ λ q随 x (本牌号 )增加而减小，呈负指数关系假设 C 实例1 产销量安排 0 y x VOR2 x=629, y=375 309.0 0 (1.3 0 ) 864.3(2.0) 飞机 x=?, y=? VOR1 x=764, y=1393 161.2 0 (0.8 0 ) VOR3 x=1571, y=259 45.1 0 (0.6 0 ) 北 DME x=155, y=987 飞机与监控台（图中坐标和测量距离的单位是 “公里 ”） 实例 2 飞机精确定位问题 飞机精确定位模型 4 2 4 2 4 )()( ),(atan2 dyyxx yyxx iii =?+? =?? θ 不考虑误差因素 超定方程组， 非线性最小二乘！ ）飞机位置坐标（要求计算： ，距离误差为记测量距离为 ，角度误差为记测量角度为 标分别为已知数据：设备位置坐 yx d i iyx ii ii , . ;3,...,1, 1,...4;),,( 44 σ σθ = = 量纲不符！ ？ 2 4 2 4 2 4 3 1 2 )()(),(atan2 ? ? ? ? ? ? ? ? ?+? + ? ? ? ? ? ? ? ? ?? ∑ = d yyxxyyxx Min i i ii x,y θ 飞机精确定位模型 44 2 4 2 444 )()( ),(atan2 σσ σθσθ +≤?+?≤? +≤??≤? dyyxxd yyxx iiiiii 考虑误差因素 Min x; Min y; Max x; Max y. 以距离为约束，优化角度误差之和（或平方和）； 或以角度为约束，优化距离误差 . 非线性规划 误差非均匀分布！ ？ ？ 仅部分考虑误差 ! 角度与距离的 “地位 ”不应不同！ 有人也可能会采用其他目标，如： 飞机精确定位模型 2 4 2 4 2 44 2 3 1 )()(),(atan2 ),( ? ? ? ? ? ? ? ? ?+?? + ? ? ? ? ? ? ? ? ??? = ∑ = σσ θ yyxxdyyxx yxEMin i i iii 误差一般服从什么分布？ 正态分布！ 不同的量纲如何处理？ 无约束非线性最小二乘模型 归一化处理！ 角度需要进行预处理， 如利用 atan2函数 , 值域 (-pi, pi) 无约束优化:最优解的分类和条件 )(xfMin x 给定一个函数 f(x),寻找 x* 使得 f(x*)最小，即 nT n xxxx ?∈= ),,,( 21 L 其中 局部最优解 全局最优解 必要条件 0),,()( 1 * ==? T xx n ffxfL x * f(x) x l xg o 充分条件 0)(,0)( *2* >?=? xfxf Hessian阵 nn ji xx f f × ? ? ? ? ? ? ? ? ?? ? =? 2 2 3 求解无约束优化的基本思路 基本思想 在 n ? 中某一点，确定一个搜索方向 及沿该方向的移动步长，得到使目标函数下降的新的点 迭代 步骤 Step 1 初始化：初始点 x 0 ，终止准则等 Step 2 迭代改进：方向 d k ，步长 α k )()( 1 kk xfxf < +kkkk dxx α+= +1 Step 3 终止检验：得到近优解或 k+1?k转 2 选择 d k ,α k 使 f 下降更快 ? 不同算法 搜索方向的选择 1 最速下降法 （梯度法） 下降方向 最速下降方向 迭代改进格式 算法特点 初始阶段改进较快，最优解附近改进较慢 kkTkkkk dxfxfdxfxf )()()()( 1 ?+=+= + 0)( <? kkT dxf )( kk xfd ??= )( 1 kkk xfxx ??= + 将 )( 1+k xf 在 k x 点作泰勒展开，只保留一阶项，有 暂不考虑搜索步 长，可设 α k =1 (负梯度方向) 2 Newton方法 特点 局部 2阶收敛 ;需计算 Hessian阵 ,它可能病态或不正定 将 f(x k+1 )在 x k 点作泰勒展开至二阶项，用 d替代 d k 牛顿方程 牛顿方向 迭代格式 )()( 2 kk xfdxf ??=? dxfddxfxfdxfxf kTkTkkk )( 2 1 )()()()( 21 ?+?+=+= + )())(( 12 kkk xfxfd ???= ? )())(( 121 kkkk xfxfxx ???= ?+ )()]([ 11 kkkk xFxFxx ?+ ′?= 比较 的牛顿法 解 0)( =xF )()( xfxF ?? 求 d使 f(x k+1 )极小 ?右端对 d导数为0 0)()( 2 =?+?? dxfxf kk 3 拟 Newton方法 目的 不计算 Hessian阵，克服病态、不正定、计 算复杂等缺陷，同时保持收敛较快的优点 回顾解方程组 F(x)=0的拟牛顿法 )()()( 11 ?? ?=? kkkkkk xFxFxxAA 满足使 kkkk AAAA 计算用迭代方法 11 ?? ?+= 思路 )()]([ 11 kkkk xFxFxx ?+ ′?= )()( 11 kkkk xFAxx ?+ ?= )(xfMin x 优化问题 )()( xFxf 相当? fGfxFxf 222 ),()( ??′? 代替阵不一定正定，构造正定相当 3 拟 Newton方法(续) kkkkkk HHHGGG ?+=?+= ++ 11 或迭代公式构造 按照 拟牛顿条件： kkkkkk fHxfxG ?=??=? ++ 11 或 )()(, 11 kkkkkk xfxffxxx ???=??=? ++ 记 )( 1 kkkk xfHxx ??= + 设在第 k步 , G k 已得到 , H k =(G k ) -1 ,可计算 )( 1112 ++++ ??= kkkk xfHxx于是有 3.1 Davidon-Fletcher-Powell(DFP)公式 kkTk kTkkk kTk Tkk k fHf HffH fx xx H ?? ?? ? ?? ?? =? )( )( )( )( 3 拟 Newton方法 (续 ) kTk TkkkkTkk kTk Tkk kTk kkTk k fx fxGGxf fx ff fx xGx G ?? ??+?? ? ?? ?? ?? ?? +=? )( )()( )( )( ) )( )( 1( 3.2 Broyden-Fletcher-Goldfarb-Shanno(BFGS)公式 kkTk kTkkk kTk Tkk k xGx GxxG xf ff G ?? ?? ? ?? ?? =? )( )( )( )( kTk TkkkkTkk kTk Tkk kTk kkTk k fx xfHHfx fx xx fx fHf H ?? ??+?? ? ?? ?? ?? ?? +=? )( )()( )( )( ) )( )( 1( 4 搜索步长的确定 线性搜索 (Line Search)确定步长方法 问题 给定 x k 和方向 d k , 确定步长 α k , 使得 )(min kk dxf α α + 优化 算法 黄金分割 (0.618)法、 Fibonacci法、 Newton 切线法、割线法、 2次或 3次插值法 等 ~一维优化问题 0=k d df α α 0)()( =+? kkkTk dxfd α 非线性最小二乘(Least Square)拟合 问题 给定 (t i , y i ), i=1,…n, 拟合一个函数 y=f(t, x)， 其中 x为待定的参数向量 , f 对 x非线性。 优化模型 )()()(min xrxrxR T x = ),()( xtfyxr iii ?= T n xrxrxr ))(),(()( 1 L= 记误差 根据目标函数是 r(x) 的二次函数的特点 构造简单算法 非线性最小二乘拟合 )()(2 xrxJR T =? SxJxJR T 2)()(2 2 +=? mmlkii xxrxr × ???=? )()( 22 讨论 ? 牛顿法要计算 Hessian矩阵，其中 S计算量大 ? 若 f 对 x线性 , 则化为线性最小二乘拟合 ,此时 S=0 特定算法考虑如何忽略或近似矩阵S。 mnji xrxJ × ??= )/()(记 )(xr 的雅各比阵为 )()()( xrxrxR T = )()( 2 1 xrxrS i n i i ?= ∑ = Gauss-Newton算法： 忽略矩阵 S )()()()( kkkkTk xrxJdxJxJ ?= f用 R代替，下降方向 d k 满足 G-N算法 收敛性依赖 f 对 x 的线性程度 , 及偏差 r的大小 非线性最小二乘拟合 )()( 2 kk xfdxf ??=?牛顿方程 )()(2 xrxJR T =? )()(2 2 xJxJR T =? Levenbery-Marquardt算法： G-N算法修正 )()())()(( kkkkkTk xrxJdIxJxJ ?=+α 其中 α k >0为修正参数 . d k 位于牛顿方向 (α k 很小 )和负梯度 方向 (α k 很大 )之间 L-M算法 非线性最小二乘拟合 )()()()( kkkkTk xrxJdxJxJ ?= 防止 J T J 出现病态 优化工具箱 基本用法： 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 RxxfMin ∈),(模型： 2 min 2 22 == + ba b y a x 其中 ：例 Exam0701.m Exam0702.m ]8,1[x )x3sinx(min 1 ∈ +：例 5 非线性最小二乘法 )()()(min xrxrxR T ubxlb = ≤≤ ),()( xtfyxr iii ?= T n xrxrxr ))(),(()( 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 xrxrxR T ubxlb = ≤≤ 3.015.247.459.3212.8914.1015.3618.1519.21c 864321.510.50.25t kt retc ? =)( 例 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 1,10,min 2 2 2 1 ==+ ba 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 222 )1()(100),(min.5 xxyyxf ?+?=例 精确解： x=y=1, f(x,y)=0 Exam0705.m 无约束优化 n x RxxfMin ∈),(模型： 采用分析梯度： GradObj='on' x=fminunc(@fun,x0,opt,…) fun.m中还要有 一般形式 function [f,g]=fun(x) )( xf? 222 )1()(100),(min.5 xxyyxf ?+?=续例 ? ? ? ? ? ? ? ? ? ???? =? )(200 )1(2)(400 2 2 xy xxyx 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 x 0 * O 的图形 222 )1()(100),( xxyyxf ?+?= exam0705plot.m Rosenbrock 函数 （香蕉函数） 算法选择 ： BFGS公式，混合 2， 3次插值，一般较好。 几个值得注意的问题 精度控制 ：对迭代次数有重大影响，应适当选择。 梯度函数 ：利用分析梯度 可能 改进算法的性能 改变初始值 由一个初值出发通常得到局部最优解， 如果函数存在多个局部最优，只有改变初值，对局 部最优进行比较，才有可能得到全局最优解。 其他算法选择 ： （详细用法请查阅 help文档） 高度非线性、不连续时可用程序 fminsearch(@fun,x0) 单变量时可用程序 fminbnd(@fun,v1,v2) 7 非线性最小二乘法 )()()(min xrxrxR T ubxlb = ≤≤ ),()( xtfyxr iii ?= T n xrxrxr ))(),(()( 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) )(xr? 算法选择 : 缺省：大规模算法（ LargeScale = 'on' ） 当 LargeScale = 'off' ： ?缺省： Levenberg-Marquardt算法； ?LevenbergMarquardt=‘off’:Gauss-Newton法 一维搜索（线搜索）步长选择与 fminunc中类似 非线性最小二乘法 )()()(min xrxrxR T ubxlb = ≤≤ exam0705lsq.m exam0705lsq_jacob_run.m 实例 1 产销量安排 ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? 280 100 22.0 1.01 2221 1211 b aa aa 已知 数据 ( ) ( ) )3020(,02.0015.0,10030 === TTT cr λ 2222221212 1112121111 )( )()(min 22 11 xcerxaxab xcerxaxabxf x x ????? ?????= ? ? λ λ ,2,1,0,,, 212211 =>??= iaabxaxabp iiiiiii ,2,1,0,,, =>+= ? icrcerq iiii x ii ii λ λ 22211121 , )()(),(max 21 xqpxqpxxz xx ?+?= 原问题 x = 23.9025 62.4977 y = 6.4135e+003 即甲产量为 23.9025，乙产量为 62.4977，最大利润为 6413.5 命令和最优解 x=fminunc(@f, x0) shili0701.m 702/502/ 22221111 ==== abxabx ，其解为 初始点 的选择 实例 1 产销量安排 ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? 22.0 1.01 2221 1211 aa aa 忽略成本及价 格中的 a 12 , a 21 2222211111 )()()(min xxabxxabxf ????= 问题简化为 作为原问题的初始点 飞机坐标 (978.31， 723.98), 误差平方和 0.6685 (<< 4) 实例2 飞机精确定位模型 2 4 2 4 2 44 2 3 1 )()(),(atan2 ),( ? ? ? ? ? ? ? ? ?+?? + ? ? ? ? ? ? ? ? ??? = ∑ = σσ θ yyxxdyyxx yxEMin i i iii 无约束非线性最小二乘模型 角度需要进行预处理， 如利用 atan2函数 , 值域 (-pi, pi) shili0702.m 其他 :非线性方程 (组 ) Rxxf ∈= ,0)(单变量方程求根： [x,fval,exitflag,output] = fzero(fun,x0,options,P1,P2,...) [x,fval,exitflag,output,jacobian] = fsolve(fun,x0,options,P1,P2, ... ) n RxxF ∈= ,0)(多变量方程组求零点： 8 2 2 2 1 0 min dCx x ? ≥ [x,resnorm,res,exitf,out,lambda]= lsqnonneg(C,d,x0,options) 数据拟合也可用 lsqcurvefit [x,resnorm,res,exitf,out,lambda,jacob] = lsqcurvefit(fun,x0,x,y,lb,ub,opt,P1,P2,…) () ∑ ?= ? ≥ i ii x ydataxdataxF ydataxdataxF 2 2 1 2 2 2 1 0 ),( ),(min 其他 :非负最小二乘 ubxlbbxAbxAts dCx x ≤≤≤≤ ? ≥ ,,.. min 2211 2 2 2 1 0 [x,resnorm,res,exitflag,output,lambda] = lsqlin(C,d,A1,b1,A2,b2,lb,ub,x0,options) 其他优化程序 Fgoalattain：多目标规划 fminimax ：极小极大 fseminf ：半无穷规划 其他优化程序 (实验 8) linprog ：线性规划 quadprog：线性规划 fmincon ：约束非线性规划 其他 :约束线性最小二乘 无约束优化实验内容 实验目的 1．掌握 Matlab 优化工具包的基本用法，对不同 算法进行初步分析、比较。 2． 练习建立和求解实际问题的无约束优化模型和 非线性最小二乘拟合模型。 内容 课上布置，或参见网络学堂