无约束最优化数学建模与数学实验实验目的实验内容
2.掌握用数学软件包求解无约束最优化问题,
1.无约束最优化基本算法,
1,无约束优化基本思想及基本算法,
4,实验作业,
3,用 MATLAB求解无约束优化问题,
2,MATLAB优化工具箱简介,
无约束最优化问题求解无约束最优化问题的的基本思想
*无约束最优化问题的基本算法返回
Rmi n nX fX?
其中 1,R Rnf?
标准形式:
求解无约束最优化问题的基本思想求解的基本思想 ( 以二元函数为例 )
1x
2x
12(,)f x x
O
1x
2x
O
5
3 1
0X
1X
2X
)( 0Xf )( 1Xf? )( 2Xf?
连续可微
Rm a x nX fX? =Rm in [ ]nX fX
多局部极小
2 9 8.0?f
0?f
2 9 8.0?f
唯一极小
(全局极小 )
221 2 1 1 2 2 1 2(,) 2 2 3f x x x x x x x x
搜索过程 2 2 2
1 2 2 1 1m i n (,) 1 0 0 ( ) ( 1 )f x x x x x
最优点 (1 1)
初始点 (-1 1)
1x 2x f
-1 1 4.00
-0.79 0.58 3.39
-0.53 0.23 2.60
-0.18 0.00 1.50
0.09 -0.03 0.98
0.37 0.11 0.47
0.59 0.33 0.20
0.80 0.63 0.05
0.95 0.90 0.003
0.99 0.99 1E-4
0.999 0.998 1E-5
0.9997 0.9998 1E-8 返回
⑴ 给定初始点 0 R nX?,允许误差
0
,令 k=0 ;
⑵ 计算
kXf?;
⑶ 检验是否满足收敛性的判别准则,
kXf
,
若满足,则停止迭代,得点 kXX?*,否则进行⑷;
⑷ 令
kk XfS
,从 kX 出发,沿 k
S
进行一维搜索,
即求
k?
,使得,
kkkkk SXfSXf

0
m i n;
⑸ 令
k
k
kk SXX 1
,k = k +1 返回⑵,
无约束优化问题的基本算法最速下降法是一种最基本的算法,它在最优化方法中占有重要地位,最速下降法的优点是工作量小,存储变量较少,初始点要求不高;缺点是收敛慢,最速下降法适用于寻优过程的前期迭代或作为间插步骤,当接近极值点时,宜选用别种收敛快的算法,
1.最速下降法(共轭梯度法)算法步骤:
2.牛顿法算法步骤:
(1) 选定初始点 0 R nX?,给定允许误差 0,令 k =0 ;
(2) 求
kXf?
, 12 kXf,检验:若
kXf
,则停止迭代,* kXX?令,否则,转 (3) ;
(3) 令
kkk XfXfS 12 ][
(牛顿方向);
(4 ) kkk SXX 1,1 kk,转回 (2),
如果 f是对称正定矩阵 A的二次函数,则用牛顿法,经过一次迭代就可达到最优点,如不是二次函数,则牛顿法不能一步达到极值点,
但由于这种函数在极值点附近和二次函数很近似,因此牛顿法的收敛速度还是很快的,
牛顿法的收敛速度虽然较快,但要求黑塞矩阵可逆,要计算二阶导数和逆矩阵,就加大了计算机的计算量和存储量,
3.拟牛顿法 为克服牛顿法的缺点,同时保持较快收敛速度的优点,利用第 k 步和第 k +1 步得到的 kX,1?kX,)( kXf?,)( 1 kXf,构造一个正定矩阵 1?kG 近似代替 )(2 kXf?,或用 1?kH 近似代替 12 ))(( kXf,将牛顿方 向改为,
1?kG 1?kS = - )( 1 kXf,1?kS = - 1?kH )( 1 kXf
从而得到下降方向,
通常采用迭代法计算 1?kG,1?kH,迭代公式 为,
BFGS ( B oryd en - F letc her - Goldfa rb - S hann o )公式
TT
1
TT
( ) ( )
( ) ( )
k k k k k k
kk
k k k k k
f f G x x G
GG
f x x G x


TT
1
TT
( ) ( )
1
( ) ( )
k k k k k
kk
k k k k
f H f x x
HH
f x f x



TT
T
( ) ( )
()
k k k k k k
kk
x f H H f x
fx


D F P ( D avid on - F letc her - Po well )公式,
TT
1
TT
( ) ( )
1
( ) ( )
k k k k k
kk
k k k k
X G X f f
GG
X f f X



TT
T
( ) ( )
()
k k k k k k
kk
f X G G X f
Xf


TT
1
TT
( ) ( )
( ) ( )
k k k k k k
kk
k k k k k
X X H f f H
HH
f X f H f



计算时可置 IH?1 (单位矩阵),对于给出的 1X 利用上面的公式进行递推,这种方法 称为 拟牛顿法,
返回
MATLAB优化工具箱简介
1.MATLAB求解优化问题的主要函数类 型 模 型 基本函数名一元函数极小 m i n F ( x ) s.t,x
1 < x < x 2
x= f mi n b nd ( ‘ F ’,x 1,x2 )
无约束极小 m i n F ( X )
X= f mi nu nc ( ‘ F ’,X 0)
X= f mi ns ea rc h ( ‘ F ’,X0 )
线性规划
m i n
Xc T
s.t,AX ≤ b
X=l in pr og (c,A,b )
二次规划
m i n
2
1
x
T
Hx + c
T
x
s.t,Ax ≤ b
X=q ua d p ro g ( H,c,A,b)
约束极小
(非线性规划)
m i n F ( X )
s.t,G ( X ) ≤ 0
X= f mi nc on ( ‘ FG ’,X0 )
达到目标问题
m i n r
s.t,F ( x ) - wr ≤ goa l
X= f go al at ta in ( ‘ F ’,x,
go a l,w)
极小极大问题
m i n max { F i ( x )}
x { F
i
( x )}
s.t,G ( x ) ≤ 0
X= f mi n i ma x( ‘ FG ’,x 0)
2.优化函数的输入变量使用优化函数或优化工具箱中其他优化函数时,输入变量见下表,
变量 描 述 调用函数
f
线性规划的目标函数 f × X 或二次规划的目标函数 X
T
× H × X + f × X 中线性项的系数向量
linprog,quadprog
fun
非线性优化的目标函数,f u n 必须为行命令对象或 M 文件、嵌入函数、或 MEX 文件的名称
fminbnd,fminsearch,fminunc,
fmincon,lsqcurvefit,lsqnonli,
fgoalattain,fminimax
H
二次规划的目标函数 X
T
× H × X + f × X 中二次项的系数矩阵
quadprog
A,b
A 矩阵和 b 向量分别为线性不等式约束:
bAX?
中的系数矩阵和右端向量
linprog,quadprog,fgoalattain,
fmincon,fminimax
Aeq,b eq
Aeq 矩阵和 beq 向量分别为 线性等式约束,
b e qXA e q
中的系数矩阵和右端向量
linprog,quadprog,fgoalattain,
fmincon,fminimax
vlb,vub X 的下限和上限向量,vlb ≤ X ≤ vub
linprog,quadp rog,fgoalattain,
fmincon,fminimax,lsqcurvefit,
lsqnonlin
X0 迭代初始点坐标 除
fminbnd 外所有优化函数
x1,x2 函数最小化的区间 fminbnd
options
优化选项参数结构,定义用于优化函数的参数所有优化函数
3.优化函数的输出变量见下表,
变量 描 述 调用函数
x
由优化函数求得的值,若 e xi t fl ag >0,则 x
为解 ; 否则,x 不是最终解,它只是迭代停 止时优化过程的值所有优化函数
fval 解 x 处的目标函数值
linpro g,qu adpr og,f goal a ttain,
fminco n,fm inim ax,l sqcu r vefit,
lsqnon lin,fmi nbnd
exitfl ag
描述退出条件,
ex i tf la g >0,表 示 目标函数收敛于解 x

ex i tf la g =0,表 示 已达到函数评价或迭代的最大次数
ex i tf la g <0,表 示 目标函 数不收敛
output
包含优化结果信息的输出结构,
It e ra ti on s,迭代次数
Al g or it hm,所采用的算法
Fu n cC ou nt,函数评价次数所有优化函数
4.控制参数选项的设置
(3) MaxIter,允许进行迭代的最大次数,取值为正整数,
选项中常用的几个参数的名称、含义、取值如下,
(1) 陈列,显示水平,取值为 'off'时,不显示输出 ;
取值为 'iter'时,显示每次迭代的信息 ;取值为
'final'时,显示最终结果,默认值为 'final'.
(2) MaxFunEvals,允许进行函数评价的最大次数,取值为正整数,
例,opts=optimset('Display','iter','TolFun',1e-8)
该语句创建一个称为选择的优化选项结构,其中显示参数设为 'iter',TolFun参数设为 1e-8.
控制参数选项可以通过函数 optimset创建或修改,命令的格式如下:
(1) options=optimset('optimfun')
创建一个含有所有参数名,并与优化函数 optimfun相关的默认值的选项结构,
( 2) options=optimset('param1',value1,'param2',value2,...)
创建一个名称为选项的优化选项参数,其中指定的参数具有指定值,所有未指定的参数取默认值,
(3) options=
optimset(oldops,'param1',value1,'param2',value2,...)
创建名称为 oldops的参数的拷贝,用指定的参数值修改 oldops
中相应的参数,
返回用 MATLAB解无约束优化问题
1,一元函数无约束优化问题,min ()fx 21 xxx
其中等式( 3)、( 4)、( 5)的右边可选用( 1)或( 2)
的等式右边,
函数 fminbnd的算法基于黄金分割法和二次插值法,它要求目标函数必须是连续函数,并可能只给出局部最优解,
常用格式如下:
( 1) x= fminbnd (fun,x1,x2)
( 2) x= fminbnd (fun,x1,x2,options)
( 3) [x,fval]= fminbnd( …)
( 4) [x,fval,exitflag]= fminbnd( …)
( 5) [x,fval,exitflag,output]= fminbnd( …)
运行结果,
xm in = 3,9 2 70 y mi n = - 0,027 9
xm ax = 0,7 854 ym ax = 0,644 8
MATLAB(wliti1)
例 1 求 x = 2 e s i nx x? 在 0 < x < 8 中的最小值与最大值,
主程序为 wliti1.m:
f='2*exp(-x).*sin(x)';
fplot(f,[0,8]); %作图语句
[xmin,ymin]=fminbnd (f,0,8)
f1='-2*exp(-x).*sin (x)';
[xmax,ymax]=fminbnd (f1,0,8)
例 2 有边长为 3m的正方形铁板,在四个角剪去相等的正方形以制成方形无盖水槽,问如何剪法使水槽的容积最大?
设剪去的正方形的边长为 x,则水槽的容积为,xx )23( 2?
建立无约束优化模型为,m in y = - xx )23( 2?,0< x <1.5
解先编写 M文件 fun0.m如下,
function f=fun0(x)
f=-(3-2*x).^2*x;
主程序为 wliti2.m:
[x,fval]=fminbnd('fun0',0,1.5);
xmax=x
fmax=-fval
运算结果为,xmax = 0.5000,fmax =2.0000.即剪掉的正方形的边长为 0.5m时水槽的容积最大,最大容积为 2m3.
MATLAB(wliti2)
命令格式为,
( 1) x= fminunc( fun,X0 );或 x=fminsearch( fun,X0 )
( 2) x= fminunc( fun,X0,options);
或 x=fminsearch( fun,X0,options)
( 3) [x,fval]= fminunc(,..);
或 [x,fval]= fminsearch(,..)
( 4) [x,fval,exitflag]= fminunc(,..);
或 [x,fval,exitflag]= fminsearch
( 5) [x,fval,exitflag,output]= fminunc(,..);
或 [x,fval,exitflag,output]= fminsearch(,..)
2.多元函数无约束优化问题标准型为,min ()FX
[3] fminunc为中型优化算法的步长一维搜索提供了两种算法,
由选项中参数 LineSearchType控制:
LineSearchType=‘quadcubic’ (缺省值 ),混合的二次和三次多项式插值;
LineSearchType='cubicpoly',三次多项式插
使用 fminunc和 fminsearch可能会得到局部最优解,
说明,?fminsearch是用单纯形法寻优,fminunc算法见以下几点说明:
[1] fminunc为无约束优化提供了大型优化和中型优化算法,
由选项中的参数 LargeScale控制:
LargeScale=‘on' (默认值 ),使用大型算法
LargeScale=‘off' (默认值 ),使用中型算法
[2] fminunc为中型优化算法的搜索方向提供了 4种算法,由选项中的参数 HessUpdate控制:
HessUpdate='bfgs'(默认值),拟牛顿法的 BFGS公式;
HessUpdate='dfp',拟牛顿法的 DFP公式;
HessUpdate='steepdesc',最速下降法例 3 min
MATLAB(wliti3)
1.编写 M文件 fun1.m:
function f = fun1 (x)
f = exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);
2.输入 M文件 wliti3.m如下,
x0 = [-1,1];
x=fminunc('fun1',x0);
y=fun1(x)
3.运行结果,
x= 0.5000 -1.0000
y = 1.3029e-10
1221 2 1 2 2( ) ( 4 2 4 2 1 ) e xf x x x x x x
1,为获得直观认识,先画出 Rosenbrock 函数的三维图形,
输入 命令,
[x,y ]= me sh gr id ( - 2,0,1,2,- 1,0,1,3 );
z= 10 0 *( y - x,^2 ),^2 + (1 - x ),^2 ;
me sh ( x,y,z)
例 4 R ose nbr o ck 函数
2 2 21,2 2 1 1( ) 1 0 0 ( ) ( 1 )f x x x x x
的最优 (极小) 解 为 x
*
= ( 1,1 ),极小值为 f
*
=0,试用不同算法(搜索方向和步长搜索)求数值最优解,
初值选为
0x
= ( - 1,2,2 ),
2,画出 R os en br oc k 函数的等高线图,输入命令,
contour(x,y,z,20)
hold on
plot( - 1.2,2,′ o ′ );
text( - 1.2,2,′ star t point ′ )
plot(1,1,′ o ′ )
text(1,1,′ solutio n ′ )
MATLAB (wliti31)
MATLAB (wliti32)
3.用 fminsearch函数求解
MATLAB (wliti41)
输入命令,
f='100*(x(2)-x(1)^2)^2+(1-x(1))^2';
[x,fval,exitflag,output]=fminsearch(f,[-1.2 2])
运行结果,
x =1.0000 1.0000
fval =1.9151e-010
exitflag = 1
output=
iterations,108
funcCount,202
algorthm,'Nelder-Mead simplex direct search '
4.用 fminunc 函数
MATLAB (wliti44)
(1)建立 M文件 fun2.m
function f=fun2(x)
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2
(2)主程序 wliti44.m
Rosenbrock函数不同算法的计算结果搜索方向步长搜索 最优解 最优值 迭代次数混合二、三次插值
( 0.999 6,0,999 2 ) 2.3109
710
1 55
B F G S
三次插值
( 1.000 1,1,000 2 ) 2.3943
810
1 32
混合二、三次插值
( 0.999 5,0,9990 ) 2.6223
710
1 51
D F P
三次插值
( 0.899 4,0,799 5 ) 0.0192 20 4
( - 1.16 34,1.36 10 ) 4.6859 20 4
( 0.944 6,0,892 0 ) 0.0031 8002
最速下降 混合 二、三次插值
( 0.99 59,0.99 16 1.8543
510
9002
单纯形法
(1.000 0,1,0000 ) 1.9151
1010
202
可以看出,最速下降法的结果最差,因为最速下降法特别不适合于从一狭长通道到达最优解的情况,
例 5 产销量的最佳安排某厂生产一种产品有甲、乙两个牌号,讨论在产销平衡的情况下如何确定各自的产量,使总利润最大,所谓产销平衡指工厂的产量等于市场上的销量,
符号说明
z ( x 1,x 2 ) 表示总利润;
p 1,q 1,x 1 分别表示甲的价格、成本、销量;
p 2,q 2,x 2 分别表示乙的价格、成本、销量;
a ij,b i,λ i,c i ( i,j = 1,2 )是待定系数,
基本假设
1.价格与销量成线性关系利润既取决于销量和价格,也依赖于产量和成本,按照市场规律,
甲的价格 p 1 会随其销量 x 1 的增长而降低,同时乙的销量 x 2 的增长也会使甲的价格有稍微的下降,可以简单地假设价格与销量成线性关系,
即,p 1 = b 1 - a 11 x 1 - a 12 x 2,b 1,a 11,a 12 > 0,且 a 11 > a 12 ;
同理,p 2 = b 2 - a 21 x 1 - a 22 x 2,b 2,a 21,a 22 > 0,且 a 22 > a 2 1,
2.成本与产量成负指数关系甲的成本随其产量的增长而降低,且有一个渐进值,可以假设为负指数关系,即,
111 1 1 1 1 1e,,,0xq r c r c
同理,
222 2 2 2 2 2e,,,0xq r c r c
模型建立若根据大量的统计数据,求出系数
b1=100,a11=1,a12=0.1,b2=280,a21=0.2,a22=2,r1=30,λ1=0.015,
c1=20,r2=100,λ2=0.02,c2=30,则问题转化为无约束优化问题:
求甲,乙两个牌号的产量 x1,x2,使总利润 z最大,
为简化模型,先忽略成本,并令 a12=0,a21=0,问题转化求
z1 = ( b1 - a11x1 ) x1 + ( b2 - a22x2 ) x2的极值,显然其解为 x1 = b1/2a11 = 50,x2 = b2/2a22 = 70,我们把它作为原问题的初始值,
总利润为,z(x1,x2)=(p1-q1)x1+(p2-q2)x2
模型求解
1.建立 M文件 fun.m,
function f = fun (x)
y1=((100-x(1)- 0.1*x(2))-(30*exp(-0.015*x(1))+20))*x(1);
y2=((280-0.2*x(1)- 2*x(2))-(100*exp(-0.02*x(2))+30))*x(2);
f=-y1-y2;
2.输入命令,
x0=[50,70];
x=fminunc('fun',x0),
z=fun (x)
3.计算结果,
x=23.9025,62.4977,z=6.4135e+003
即甲的产量为 23.9025,乙的产量为 62.4977,最大利润为 6413.5.
MATLAB(wliti5)
返回实验作业
1,求下列函数的极小点,
1)
21232221 18294 xxxxxXf;
2)
21212221 22
2
3
xxxxxxXf;
3)
2241 21 xXf
,
第 1 ),2 )题的初始点可任意选取,
第 3 )题的初始点取为
T0 0,1X?
,
2,梯子长度问题一楼房的后面是一个很大的花园,在花园中紧靠着楼房有一个温室,高 3 m,温室伸入花园 2 m,
温室正上方是楼房的窗台,清洁工打扫窗台周围,
他得用梯子越过温室,一头放在花园中,一头靠在楼房的墙上,因为温室是不能承受
a 梯子压力的,所以梯子太短是不行的,
现清洁工只有一架 7 m 长的梯子,
b 你认为它能达到要求吗? 能满足要求的梯子的最小长度为多少?
3,陈酒出售的最佳时机问题某酒厂有批新酿的好酒,如果现在就出售,
可得总收入
0R
=50 万元 ( 人民币 ),如果窖藏起来待来日 ( 第 n 年 ) 按陈酒价格出售,第 n 年末可得总收入
6
0 e
nRR?
( 万元 ),而银行利率为 r =0,0 5,试分析这批好酒窖藏多少年后出售可使总收入的现值最大,( 假设现有资金 X 万元,将其存入银行,到第 n 年时增值为
()Rn
万元,则称 X 为
()Rn
的现值,) 并填下表,
第一种方案:将酒现在出售,所获 50 万元本金存入银行;
第二种方案:将酒窖藏起来,待第
n
年出售,
( 1 ) 计算 15 年内采用两种方案,50 万元增值的数目并填入表 1,2 中;
( 2 ) 计算 15 年内陈酒出售后总收入
()Rn
的现值填入表 3 中,
表 1 第一种方案第 1 年 第 2 年 第 3 年 第 4 年 第 5 年第 6 年 第 7 年 第 8 年 第 9 年 第 10 年第 11 年 第 12 年 第 13 年 第 14 年 第 15 年表 2 第二种方案第 1 年 第 2 年 第 3 年 第 4 年 第 5 年第 6 年 第 7 年 第 8 年 第 9 年 第 10 年第 11 年 第 12 年 第 13 年 第 14 年 第 15 年表 3 陈酒出售后的现值第 1 年 第 2 年 第 3 年 第 4 年 第 5 年第 6 年 第 7 年 第 8 年 第 9 年 第 10 年第 11 年 第 12 年 第 13 年 第 14 年 第 15 年