优化设计方案

工学院课程设计报告

课程名称：优化计算方法课程设计

指导老师：

班级：交运174班

姓名：

学号：

学期：20 18 —20 19 学年第 2 学期

南京农业大学工学院教务处印

目录

1、抛物线法与黄金分割法的比较 (2)

1.1黄金分割法原理 (2)

1.2抛物线法原理 (3)

1.3黄金分割与抛物线法求解问题 (4)

1.4使用matlab求解过程 (4)

1.5编写与运行制图程序 (6)

2、matlab优化工具箱的使用 (7)

2.1求解一个约束非线性问题 (7)

2.2运行优化 (9)

2.2.1使用优化应用程序最小化Rosenbrock的功能 (9)

2.3最小化Rosenbrock在命令行中的功能： (10)

3、Matlab工具箱——Fmincon函数的使用 (10)

3.1算法简介 (10)

3.2．算法相关句法 (11)

3.3例子——非线性约束 (12)

3.4利用fmincon函数解决实际问题 (14)

3.4.1编程过程 (14)

3.5结果分析 (15)

3.5.1运行行优化 (15)

优化计算方法课程设计

交通运输专业 曹福 指导老师 朱磊

1、抛物线法与黄金分割法的比较

1.1黄金分割法原理

黄金分割法也称为0.618法，其基本思想是：通过试探点函数值的比较，使包含极小点的搜索区间不断缩小，该方法仅需要计算函数值，使用范围广，使用方便。

设

)()(k k sd x f s +=Φ

式中：)(s Φ为搜索区间[0a ，0b ]上的单峰函数。

在第i 次迭代时，搜索区间为[i a ，i b ]，取两个试探点为i p ，i q ∈[i a ，i b ]且i p )(i q Φ，则令i i p a =+1，i i b b =+1。

要求两个试探点i p 和i q 满足下述两个条件： (1) [i a ，i q ]与[i p ，i b ]的长度想同，即i b -i p =i q -i a ； (2) 区间长度的缩短率相同，即1+i b -1+i a =t (i b -i a )。 从而得

i p =i a +(1-t )(i b -i a )，i q =i a +t (i b -i a )

考虑情况(1)，此时新的搜索区间为[1+i a ,1+i b ]=[i a ,i b ]。 为了进一步缩短搜索区间，取新的试探点1+i p 和1+i q ，得到

1+i q =1+i a +2t (i b -i a )。

若令

2t =1-t ，t >0,

则

1+i q =i a +(1-t )(i b -i a )=i p

这样新的试探点1+i q 就不用重新计算。类似的，对于情况(2)，也有相同的结论。 解方程2t =1-t ，t >0得到区间长度缩短率为

t =

0.6182

1

-5≈ 1.2抛物线法原理

抛物线法也叫二次插值法，其基本思想是：在搜索区间中不断地使用二次多项式去近似目标函数，并逐步用插值多项式的极小点去逼近线搜索问题

)()(min 0

k k s sd x f s +=Φ>

的极小点。

设已知三点

0s ，h s s +=01，h s s 202+= ()0>h

处的函数值0Φ，1Φ，2Φ，且满足

01Φ<Φ，21Φ<Φ

上述条件保证了函数在区间[0s ,2s ]上是单峰函数。则满足上述条件的二次拉格朗日插值多项式为

22

101

22002212))(())((2)

)(()(Φ--+Φ---Φ--=

h s s s s h s s s s h s s s s s q )(s q 的一阶导数为

22

011220022

122222)('Φ--+Φ---Φ--=

h

s s s h s s s h s s s s q 令)('s q =0，解得

h

s

h

s

h

s

h

s

h

s

s

s

s

s

s

s

s

+

=

Φ

+

Φ

-

Φ

Φ

+

Φ

-

Φ

+

=

Φ

+

Φ

-

Φ

Φ

+

+

Φ

+

-

Φ

+

=

Φ

+

Φ

-

Φ

Φ

+

+

Φ

+

-

Φ

+

=

2

1

2

1

2

1

2

1

2

1

2

1

1

2

2

1

)

2

(2

)

4

3(

)

2

(2

)

2(

)

2

2(2

)

3

2(

)

2

(2

)

(

)

(2

)

(

式中

)

2

(2

)

4

3(

2

1

2

1

0>

Φ

+

Φ

-

Φ

Φ

+

Φ

-

Φ

=

h

h

又因为)

(s

q的二阶导数为

2

2

)

(''

2

2

1

2

2

2

1

2

0>

Φ

+

Φ

-

Φ

=

Φ

+

Φ

-

Φ

=

h

h

h

h

s

q

故)

(s

q为凸二次函数，从而

m in

s是)

(s

q的全局极小点。

1.3黄金分割与抛物线法求解问题

利用抛物线法和黄金分割法求函数

()x

x

x tan

2

32-

=

Φ

在[0，1]上的极小点。取容许误差-4

10

=

ε，-5

10

=

δ。

1.4使用matlab求解过程

步骤一：找到目标函数、梯度、黄金分割和抛物线的.m文件（如图1）

图1 目标函数、梯度、黄金分割和抛物线程序

步骤二：导入目标函数、梯度、黄金分割和抛物线程序（如图2）

图2

步骤三：输入指令运行所得结果

phi=@(x)3*x^2-2*tan(x);

[i,s,phis,ds,dphi,G]=golds(phi,0,1,1e-4,1e-5)%黄金分割法；

i =

21

s =

0.3895

phis =

-0.3658

ds =

6.6107e-05

dphi =

2.7880e-09

G =

0 0.3820 0.6180 1.0000

0 0.2361 0.3820 0.6180

0.2361 0.3820 0.4721 0.6180

0.2361 0.3262 0.3820 0.4721

0.3262 0.3820 0.4164 0.4721

0.3262 0.3607 0.3820 0.4164

0.3607 0.3820 0.3951 0.4164

0.3820 0.3951 0.4033 0.4164

0.3820 0.3901 0.3951 0.4033

0.3820 0.3870 0.3901 0.3951

0.3870 0.3901 0.3920 0.3951

0.3870 0.3889 0.3901 0.3920

0.3870 0.3882 0.3889 0.3901

0.3882 0.3889 0.3894 0.3901

0.3889 0.3894 0.3896 0.3901

0.3889 0.3892 0.3894 0.3896

0.3892 0.3894 0.3895 0.3896

0.3894 0.3895 0.3895 0.3896

0.3894 0.3894 0.3895 0.3895

0.3894 0.3895 0.3895 0.3895

0.3895 0.3895 0.3895 0.3895

>> phi=@(x)3*x^2-2*tan(x);

[i,s,phis,ds,dphi,S]=qmin(phi,0,1,1e-4,1e-5)%抛物线法；

i =

5

s =

0.3895

phis =

-0.3658

ds =

5.9671e-05

dphi =

7.2670e-09

S =

0 0.5069 0.3938 0.3895 0.3895

步骤四：分析结果

根据以上运算结果得到黄金分割法和抛物线法求单变量函数极小点的数值结果的对比如下

方法参数迭代次数（i）极小点（*s）

极小值

（()*

Φs）i

i

a

b-的值

()()

i

i

a

bΦ

-

Φ

的值

黄金分割法210.3895-0.36586.6107 ×10-52.7880 ×10-9抛物线法 5 0.3896 -0.3685 5.9671×10-5 7.2670×10-9

表1结果对比

1.5编写与运行制图程序

golds函数中：

function [i,s,phis,ds,dphi,G,Y]=golds(phi,a,b,epsilon,delta)

%在golds函数的输出列表中加入矩阵Y。

if(phip<=phiq)

b=q; phib=phiq; q=p; phiq=phip;

h=b-a; p=a+(1-t)*h; phip=feval(phi,p);Y(i,:)=[i,phip];

else

a=p; phia=phip; p=q; phip=phiq;

h=b-a; q=a+t*h; phiq=feval(phi,q);Y(i,:)=[i,phiq];

end

%在循环体中，把第i次迭代的最优解以及i本身的值分别赋给矩阵Y的第2和第1列。ds=abs(b-a); dphi=abs(phib-phia);plot(Y(:,1),Y(:,2),'--');hold on;

%以矩阵Y第1列元素为自变量，第2列元素为因变量，线型为“--”，绘制函数图像，并保持函数图像持续显示。

qmin函数中：在qmin函数输出列表中也加入矩阵Y。

bars=s0+barh;barphi=feval(phi,bars);Y(i,:)=[i,barphi];plot(Y(:,1),Y(:,2));lege nd('golds','qmin');

%在循环体中，把第i次迭代的最优解以及i本身的值分别赋给矩阵Y的第2和第1列。以矩阵Y第1列元素为自变量，第2列元素为因变量，绘制函数图像。并添加图例以示区别。

结果分析：

每次迭代得到的目标函数值可以画成折线图表示，用不同的方法求解问题得到的折线图也不一样，下图将两种方法得到的折线合并到一张图中（如图2）：

图2

抛物线法收敛速度快，总的迭代次数为5次，第1次次迭代目标函数值下降的快，之后的迭代目标函数值下降速度缓慢。

黄金分割法收敛速度慢，迭代次数为21次，前4次迭代的目标函数值波动，后基本接近极小值，并维持平衡。

2、matlab优化工具箱的使用

2.1求解一个约束非线性问题

典型优化问题

此示例显示了如何使用优化工具箱求解器解决受约束的非线性问题。该示例演示了典

型的工作流程:创建目标函数、创建约束、解决问题和检查结果。

问题表述:罗森布鲁克的功能

考虑最小化罗森布鲁克函数的问题

f(x)=100(x(2)?x（2）^1)^2+(1?x(1))2,

在单位圆盘上，即半径为1的圆盘以原点为中心。换句话说，找到x，使函数f(x)在x(1)^2+x(2)^2≤1的集合上最小化。这个问题是具有非线性约束的非线性函数的最小化。

注意：

罗森布鲁克函数是最优化中的标准测试函数。它在[1，1]点处具有唯一的最小值0。对于一些算法来说，找到最小值是一个挑战，因为函数在一个深深弯曲的山谷中有一个浅的最小值。这个问题的解决方案不在点[1，1]，因为该点不满足约束条件。

该图显示了单位圆盘中罗森布鲁克函数的两个视图。垂直轴是对数刻度的；换句话说，该图显示了log(1+f(x))。等高线位于曲面图下方。

图3

生成图形的代码

函数f(x)称为目标函数。目标函数是我们想要最小化的函数。不等式x(1)^2+x(2)^2≤1称为约束。约束限制求解器搜索最小值的x集合。您可以有任意数量的约束，即不等式或等式。

所有优化工具箱优化函数都将目标函数最小化。要最大化函数f，应用优化例程来最小化–f。

在命令行中，输入（如图4）：

图4

在MATLAB编辑器中，输入（图5）：

图5

使用名称保存文件rosenbrock.m。

创建约束条件：

创建一个unitdisk.m包含以下代码的文件（图6）：

图6

保存文件unitdisk.m。

2.1运行优化

使用优化应用程序最小化Rosenbrock的功能：

通过optimtool在命令行输入来启动优化应用程序。

在Objective function框中，输入@rosenbrock。@字符表示文件的函数句柄（MATLAB）rosenbrock.m。

在Start point框中，输入[0 0]以指定fmincon开始搜索最小值的初始点。

在Nonlinear constraint function框中，输入@unitdisk，函数句柄unitdisk.m。

在Options窗格的Display to command window窗口（位于窗格底部）下，iterative 从Level of display列表中进行选择。（如果没有看到该选项，请单击Display to command window）

在Problem Setup and Results窗格的Run solver和view results下，单击Start。

2.2最小化Rosenbrock在命令行中的功能：

创建选择迭代显示和interior-point算法的选项。

运行fmincon具有options结构的求解器，报告x最小化器的位置和fval目标函数获得的值。

3、Matlab工具箱——Fmincon函数的使用

3.1算法简介

利用Fmincon函数进行非线性规划求解。

查找由指定的问题的最小值：

变量定义解释：

b和beq是向量，a和aeq是矩阵，c（x）和ceq（x）是返回向量的函数，f（x）是返回标量的函数。f（x）、c（x）和ceq（x）可以是非线性函数。x、lb和ub可以作为向量或矩阵传递；请参见矩阵参数。

3.2．算法相关句法

其中句法说明：

1.x=fmincon（fun，x0，a，b）从x0开始，并试图找到fun中描述的函数的极小x，服从线性不等式a*x≤b。x0可以是标量、向量或矩阵。

2.x=fmincon（fun，x0，a，b，aeq，beq）根据线性等式aeq*x=beq和a*x≤b将fun最小化。如果不存在不等式，则设置a=[]和b=[]。

3.x=fminco（fun，x0，a，b，aeq，beq，lb，ub）定义了x中设计变量的一组下限和上限，因此解决方案始终在lb≤x≤ub范围内。如果不存在相等项，则设置AEQ=[]和BEQ=[]。如果x（i）在下面未绑定，则设置lb（i）=-inf；如果x（i）在上面未绑定，则设置ub（i）=inf。

4.x=fmincon（fun，x0，a，b，aeq，beq，lb,ub,nonlcon）subjects the minization to the nonlcon nonlinear in equalities C（X）or equalities CEQ（X）defined in nolcon.FMINCON Optimizes such that C（X）≤0 and CEQ（X）=0.如果没有边界存在，则设置LB=[]或UB=[]。

5.x=fmincon（fun，x0，a，b，aeq，beq，lb,ub,nonlcon，options）minimizes with the optimization options specified in options.使用选择性设置这些选项。如果没有非线性不平等或平等限制，则设置nolcon=[]。

6.x=fmincon（Problem）发现问题的最小值，其中问题是输入贡献中描述的结构。通过从优化应用中导出问题来创建问题结构，如在导出工作中描述问题。

7.[x,fval] = fmincon(___),对于任何语法，返回解决方案x中目标函数fun的值。

8.[x,fval,exitflag,output] = fmincon(___) 此外，还返回一个值exitflag，该值描述fmincon的退出条件，并返回一个结构输出，其中包含有关优化过程的信息。

9.[x,fval,exitflag,output,lambda,grad,hessian] = fmincon(___)另外还返回：

lambda-包含解决方案x处拉格朗日乘数的字段的结构。

Grad-解决方案x的函数梯度。

Hessian-解决方案x的函数所在。参见fmincon hessian。

3.3例子——非线性约束

求受非线性约束的函数的最小值。找出罗森博克函数在一个圆内最小化的点，也受约束。

1.函数问题：fun = @(x)100*(x(2)-x(1)^2)^2 + (1-x(1))^2;

2.在, 区域内查看

转化成语言为：

lb = [0,0.2];

ub = [0.5,0.8];

3.也可以在半径为1/3的[1/3,1/3]处居中的圆内查看。将以下代码复制到名为circlecon.m的matlab?路径上的文件中。

转化成语言为：

% Copyright 2015 The MathWorks, Inc.

function [c,ceq] = circlecon(x)

c = (x(1)-1/3)^2 + (x(2)-1/3)^2 - (1/3)^2;

ceq = [];

4.由于没有线性约束，因此将这些参数设置为[]。

转化成语言为：

A = [];

b = [];

Aeq = [];

beq = [];

5.选择满足所有约束的初始点：x0 = [1/4,1/4];

6.在命令行窗口打上如下语言，得出结果问题解决

nonlcon = @circlecon;

x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)

7.找到满足约束的局部最小值。由于目标函数在可行方向，在最优性公差的默认值内，并且在约束公差的默认值内满足约束。

结果为x = 0.5000 0.2500

3.4利用fmincon 函数解决实际问题

问题：

122

12122121212min ()(42421)

1.510

x f x e x x x x x x x x x x x =++++--≤-??

≥-?目标函数约束条件

3.4.1编程过程

（1）新建M 文件函数fun1，定义函数f=fun1（x ），代入目标函数：

f=exp(1)^x(1)*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1)，以fun1命名M 文件并进行保存。

程序如下所示：

function f=fun1(x)

f=exp(1)^x(1)*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1); end

（2）新建M 文件函数fun2，定义函数[C,Cep]=fun2(x)，由题可知存在两个非线性不等式，所以分别用C(1)、C(2)表示对应非线性不等式；由于不存在非线性等式所以用[]表示参数。将题中约束条件更改为C （x ）≤ 0的格式： x(1)*x(2)-x(1)-x(2)+1.5≤0； -x(1)*x(2)-10≤0并编入文件中，以fun2命名M 文件并进行保存。

程序如下所示：

function [C,Cep]=fun2(x)

C(1)=x(1)*x(2)-x(1)-x(2)+1.5;% 非线性约束，若不止1个，则C(1),C(2),... C(2)= -x(1)*x(2)-10; % 非线性约束，若不止1个，则用C(1),C(2),... Cep=[];% 等式约束, 若没有等式约束，可令C=[]; end （3）确定满足约束条件的初始值x0=[-1,2],由于不存在线性不等式、线性等式、非线性等式以及x(1)和x(2)的边界，因此A 、b 、Aeq 、beq 、lb 、ub 参数均用[]表示。nonlcon 部分存在于函数fun2中。将上述变量导入命令行窗口 程序如下所示： x0 = [-1,2]; A = []; b = []; Aeq = [];

beq = [];

lb = [];

ub = [];

[x,fval] = fmincon('fun1',x0,[],[],[],[],[],[],'fun2')

（4）得出计算结果为：

x =

-9.5474 1.0474

fval =

0.0236

3.5．结果分析

X(1)= -9.5474 X(2)=1.0474

题目中约束条件为：

x(1)*x(2)-x(1)-x(2)≤-1.5；-10≤x(1)*x(2)

代入数值可得：-1.49994676≤-1.5,-10≤9.99994676；可知满足约束条件。

代入目标函数得：min f(x)=0.0236

4、实验总结与收获

通过本次实验课程设计我学习到了MATLAB软件的简单使用方法，对MATLAB软件有了一个基本的了解，更进一步更清晰的认识了黄金分割法和抛物线法两种方法去求单变量函数极小点的值的区别，同时在这个过程中我的自我探索能力和自学能力得到了一定程度的提升。

南京农业大学工学院课程考核情况表

（实习、课程设计）

20 —20 学年第学期课程名称：

班级：

姓名：

学号：

学生表现：

出勤

学习态度

学习主动性

总体表现

报告（说明书）评语：

优点

不足

总体评价

（优良中差）

成绩：

教师签名：

日期：

（此页由老师填写）