优化设计课后题和大作业

机械优化设计作业

第一、机械优化设计课后练习

第一章机械优化的基本概念与数学模型

1-1.优化设计问题的数学模型是由哪几部分组成的？其一般表达形式是什么？

答：由三部分组成：设计变量、约束条件和目标函数。

优化设计一般表达形式是：

山]屋

Findx=[A-],X2,—RN-—设计变量

minf(x)--------目标函数

sddel

s.t.队X)=00=1〜4)--------约束条件

g」(x)wo0=小1〜J)—约束条件

式中：c={xlXL<X<XU}

1-2.建立优化设计问题数学模型的一般步骤及其需要注意的问题是什么？

答：建立优化设计问题数学模型的一般步骤为：

(1)选取设计变量

(2)建立目标函数

(3)确定约束条件

其注意事项：

(1)设计变量

在设计过程中选择的设计变量必须都是独立变量，有明显依赖关系；设计变量的选取与优化层次

及优化问题的提法有关；设计变量的数目要适当；设计变量有显著且能直接调整控制参数。

(2)约束条件

周密分析、合理确定约束条件，从客观实际出发，且能表为设计变量的约束函数的限制确定为约

束；各约束条件应当是独立而不矛盾；要特别注意那些对优化效果确有影响，确有限制作用的约束，应

注意它们是否可以适当放松以达更好优化效果。

(3)目标函数

目标函数可能是多种，具体选哪个取决于对设计的具体要求和客观条件；根据工程实际选定最重

要的为优化目标；考虑当前设计方案的实际情况；同时应考虑该指标是否容易给出数学表达式，常常以

多目标优化使用更符合实际。

1-3.优化设计问题的求解方法有哪几类？迭代法的基本思想及特点是什么？

答：①优化设计问题的求解方法分为两大类：简单优化问题的求解和数值迭代法。

(1)简单优化问题的求解方法：

a、解析法：适用于形式简单、容易求导，可直接写出数学模型显式表达式的、不带或仅带简单等

式约束的优化问题，可通过高等数学的极值条件解方程求解。

b、图解法：NW2维情况，通过作图求解，简单直观。

(2)数值迭代法——现代优化设计的基本方法

a、数学规划法：根据函数及其导数的局部性态决定迭代方向和步长。

迭代通式:X（z）=X㈤+泗）5（。

b、准则法：多用于结构优化一复杂结构优化（隐函数）。

迭代通式为X（*+D=°（X%（未必收敛）

②迭代法的基本思想：根据目标函数的变化规律，以适当的步长沿着能使目标函数下降的方向，

“步步逼近”，“步步下降”或“步步登高”，逐步向目标函数值的最优点进行探索，逐步逼近或直至达到目

标函数的最优点。

③迭代法的特点：具有简单的逻辑结构并能进行反复的同样的算术计算；最后得出的是逼近精确解

的近似解。

1-4.欲造容积为V的长方形无盖水箱，如何选定其长、宽、高，才能使用料最少，写出数

学模型。

/

解：选取长⑴，宽(d),高(h)为设计变量：X=[l,d,h^=[xl,x2,x3]

设计目标：用料最省

f(x)-S-2x(h-d+h-l)+d-I—>min

约束条件

h(X)-xtx2x3-V

g|(X)=-x,<0

g2(X)=-*2-0

g3(X)=-x3<0

数学模型为：

FindX=r€R3

min/(X)=X|X2+2(x2x3+%,%3)

h(X)=x1x2x3=V

gi(X)=*0

Stg2(X)=f40

g3(X)=-x3<0

第二章优化设计的最优性条件

2-1.梯度与海森阵的表达与意义是什么？梯度与方向导有何关系？

答：（1）梯度的表达：为一个N维列阵

W、（x）=［@更…里］

［出dx2dxNJ

意义：它是目标函数/（X）对各个设计变量的一阶偏导数所组成的列向量，方向就是函数变化率

2

最大的方向，用来研究多元函数沿各坐标方向的变化率。

（2）海森阵的表达：为二阶导数

d2fd2fd2f

dxidx2dxxdxN

d2fd2fd2f

▽2/(X)=〃(X)=,2►

dx2dxl8X2dx2dxN

57d2fd2f

dxdxdxdxdx2

NxN2NNxN

海森阵是多元函数关于诸设计变量的二阶导数矩阵，是对称方阵。

（3）梯度与方向导数的关系：多元函数在某方向上的方向导数是梯度在该方向上的投影。

2-2、求几种特殊函数的梯度与海森阵：

线性函数：B'X+c；二次型函数：XZX；一般二次函数：X'AX+B'X+c

解：（1）、线性函数：BTX+C

设X=卜|,%2,…,xj，也,…,九］则

•X

T2

/(X)=BX+c={b],b2,--,bN].+c=b}x}+b2x2+--+bNxN+c

_XN_

所以且31=a

oxi

rr

故：Vf(x)=V(BX+C)=[bI,b2,-,^]=B

因为誓

0,所以海森阵H（X）=0

（2）、二次型函数：X’AX

a\\a\2a\N

、几a2\a22%为对称阵

设：AA=

_aN\aN2aNN_

a\\a\2a\N

N

a

T2NLaxx

f(X)=XAX=[Xlx2…xj的陶X,v]'=Y'i：i

51

aN\aN2aNN

3

0(X)=2Z。卢=2AX

则海森阵为：H(X)=2A

(3)一般二次函数：XTAX+BTX+C

由上题结果，设A为对称阵，

则梯度为：Y/(X)=2AX+8,

海森阵为：〃(X)=2A

2-3、多元函数的无约束极值、等式约束极值及不等式约束极值的必要条件的具体形式是

什么？充分条件是什么？

答：(一)多元函数的无约束极值的必要条件是：设多元函数/(X)在X*处有一阶及二阶连续偏导数

(1)必要条件：V/(X*)=O

(2)充分条件：

①V/(X")=O

OH(X*)为正定或者是负定，即二阶导海森阵正定或者负定，即对任何非零N维向量丫有

yr”(x*)y〉o或＜0,正定时有极小值，负定时有极大值。

(二)多元函数的等式约束极值必要条件为Largrange条件:

在极值点X*,满足1),2),3)三个条件：

生=笠+£dh.

1).少项二°(n=l~N)

私“dxn

2).%(X)=0。=0〜J)

3).正则条件；要求N4(X*),…,V/(X*)］秩为J。

充分条件：在Largrange条件基础上，满足目标函数的高阶条件与凸规划等理论判别。

(三)多元函数的不等式约束极值必要条件：

1).W^VAX与+Nzvg/X*^。

>1

2).>00=1〜J)

3).4g/X*)=0(j=T〜J)

4).gj(X*)40。=1〜J)

另外还要满足正则条件，即在X*点诸临界约束的梯度线形无关，相互独立。

充分条件：

当f(X)为凸函数、可行域D为凸集的凸规划问题时、必要条件也就是他的充分条件，因此充分

条件为:

4

f(X)为凸函数、可行域D为凸集

▽L=W(X*)+Z%Vgj(X*)=O

4NO)=(1〜J)

48(X*)=0，=。〜J)

5).g/X*)40j=(l-J)

另外还要满足正则条件，即在X*点诸临界约束的梯度线形无关，相互独立。

2-4、求/(X)=4+4.5%]-4%2+X；+2X22-2%1%2+%4-。―2的极值点及其性质，图解验

证。

解:由W(X*)=0得:

4.5+2x]-2X2+4X「-4X]X2=0

dx1

df

—4+4%2—2X1—2工]—0

得三个解及相应的函数值，如下：

X1*=[1.9413.854『fx2*=[-1.0531.028/J%3*=[0.61171.4929y

/(%1*)=0.9855'(/(X2*)=-0.5134'=2.8091

求导得海森阵：

d2f的/

dx：dxdx2+12.V1~—4%2—2—4X1

H=V2/(X)t2

d2fd2f-2-4x,4

dx2dx}dr；

317938-9764

将X1*=[1.9413.854丫代入上式得：H='是正定阵，因此

1L」|_-9.7644_

X,*=[1.9413.8541是局部极小值点

「1119372212-

将X,*=[—1.0531.028『代入上式得：H=■也是正定阵，因此

\_2.2124

7

X2*=[-1.0531.028]也是局部极小值点。

r4「0.51896-4A469

WX*=[0.61171.49291代入上式得：H=既非正定阵，也非负定阵，

3-4.44694

5

因此乂3*=［0.61171.4929『是鞍点。

又由于〃X1*)>/(X2*),所以X2*=［-L0531.0281是最小点。

在MATLAB中画出如下图形：

从图上看出［0.6117L4929F是鞍点，［-1.0531.028『和［1.9413.854『是极小点，然后他

们的函数值，可以看出［-1.0531.028『是最小点。与计算结果相符合。

r

X=[x1,x2,x3]

22

minf(X)=x"+x2+x3

gl(X)=5-xi-x2-x3<0

g2(X)=2-x2x3<0

g3(X)=F<0

g4(X)=-x2<0

g5(X)=2-x3<0

检查(2,1,2),(4/3,2/3,3),(3/2,3/2,2)三点的Kuhn-Tucker条件。

解：不等式约束极值Kuhn-Tucker条件为:

(1)、VL=0(X*)+Z乙Vgj(X*)=0

;=1

(2)、>0j=(l〜J)

6

(3)、%gj(X*)=Oj=(l〜J)

(4)、g/X*)40，=（1〜J）

另外还要满足正则条件,即在X*点诸临界约束的梯度线形无关，相互独立。

先判定是否为凸函数

22+2

/(X)=X1+X2X3

2222

8f/dxt=2xi,8f/dX|=2,df/dx\dx2=0,df/dX\dx3=0

df/dX2=2x2，62f/dX22=2,82//5X2SXI=0>d2f/dx?。X3=0

2222

8f/dX3=2x3，Sf/dx3=2,df/dx3dxi=O,8f/dx3Sx2=0

经计算H(x)的各阶顺序主子式大于等于0,BU：H(x)半正定。

；•/(x)为凸函数。

同理可知：g(x)也为凸函数。

再求出f(X)和gj(X)在X*处的梯度：

Vf(X*)=[2x,2X22XJ

▽gi(X*)=[-l-1-l]r

▽g2(x*)=[o-x3-X2Y

▽g3(X*)=[-l00]r

▽g4(X*)=[0-10]r

▽g5(X*)=[00-l]r

①将X1*=(2,1,2)代入

gi(X*)=0,g2(X*)=0,g3(X*)=-2,g4(X*)=-l,g5(X*)=0

gi(X*)、g2(X*)和gs(X*)为起作用的约束。

TTT

贝IJVgl(X*)=[-l,-1,-1],Vgl(X*)=[0,-2,-1]Vg5(X*)=[0,0,-1]

V/(X*)=[4,2,4]T

J

代入VL=Yf(X*)+E/Vgj(X*)=00=1、2、5)即：

J=1

4-100

+%5

2+A1-1+42-200

4-1—1—1

解得：4=4,4=—i，4=i

7

因为友<0,不满足aN0,所以对于点(2,1,2)不存在满足Kuhn-Tucker条件的讥②将X2*=(4/3,

2/3,3)代入

gi(X*)=0,g2(X*)=0,g3(X*)=-4/3,g4(X*)=-2/3,g5(X*)=-l

/.gl(X*),g2(X*)为起作用的约束。

TT

Vgl(x*)=[-1,-1,-1]，Vg2(X*)=[0,-3,-2/3]▽/(X*尸[8/3,4/3,6]T代入

VL=Vf(X*)+Z4Vgj(X*)=0即:

无解，所以对于点(4/3,2/3,3)不存在满足Kuhn-Tucker条件的4。

③将X2*(3/2,3/2,2)代入

gi(X*)=O,g2(X*)=-l,g3(X*)=-3/2,g4(X*)=-3/2,g5(X*)=0

g|(X*)与g5(X*)为起作用的约束。

TTT

Vgl(X*)=[-l,-1,-1],Vg5(X*)=[0,0,-1],V/(X*)=[3,3,4]

代入VL=V/(X*)+f/l尸g/X*)=0(j=l、5)叫

A，i=3>0入5=l>0

当入=[3000时点(3/2,3/2,2)满足K—T条件。

第三章无约束规划的解法

3-1、用黄金分割法求目标函数f=x(x+2)的最优解，初始区间为［-3,5］,误差£不大于0.05。

解：(1)a=-3,b=5

N>M［”二川=10.54545，得11次迭代能得到理想解。

In0.618

="+0.382(匕一a)=0.056

在［-3,5］内取点

=a+0.6183-a)=1.944

8

则

/(x,(l))=0.115136</(x,<2))=7.667136

所以令b==1.944,即新的区间为[-3,1.944]

(2)0(1)

(2)X2=x/=0.056x2=a+0.382(&-a)=-1.1114

则

(2)(，)

/(X2)=0.115136>/(X2)=-0.9876

所以令6⑵=0.056,即新的区间为[-3,0.056]

(3)9⑵=4")=一111149⑴=a+0.382(/>—a)=—1.833

则

(1)

/(当⑵)=-0.9876</(x3)=-0.3061

所以令a=X°)=-1.833,即新的区间为[-1.833,0.056]

(4)、匕⑴=与⑵=一】1114乙⑵=a+0.618(b—a)=—0.6656

则

<2)(1)

/(X4)=-0.8882>/(x4)=-0.9876

所以令/?=/⑵=—0.6656,即新的区间为[-1.833,-0.6656]

(2)0)

(5)、x5=x4=-1.1114X5⑴=a+0.382(6-。)=-1.3871

则

(l)(2)

y(x5)=-0.8502>y(x5)=-0.9876

所以令a=/⑴=—13871,即新的区间为[-1.3871,-0.6656]

(1><2)

(6)、x6=X5=-1.1114*6⑵=a+0.618(/?—a)=—0.9412

则

(2)m

/(X6)=-0.9965<f(x6)=-0.9876

所以令a=4⑴=-1.1114,即新的区间为[-1.1114,-0.6656]

(1)<2)<2)

(7)、x7=x6=-0.9412x7=a+0.618(Z?-a)=-0.8359

则

(l)(2)

/(x7)=-0.9965</(X7)=-0.9731

9

所以令b=X7⑵=—0.8359,即新的区间为[-1.1H4,-0.8359]

(8)、/⑵=七⑴=—0.94124⑴=a+0.382(A—a)=—1.0062

则

<2)

/(X8)=-0.9965>/(%/'))=-1.0000

所以令6=4⑵=—0.9412,即新的区间为[-1.1114,-0.9412]

(2)(l)

(9)、x9=xg=-1.0062X9⑴=。+0.3823-“)=-1.0464

则

(2)(l)

/(X9)=-1.0000</(x9)=-0.9978

所以令a=/⑴=-1.0464,即新的区间为[-1.0464,-0.9412]

(10),巧0⑴=/⑵=-10062玉0⑵=a+0.618(b-a)=—0.9814

则

/(均⑵)=-0.9997>/(小⑴)=-1.0000

所以令匕=占0⑵=一0.9814,即新的区间为[-1.0062,-0.9814]

(11)、x/2)=演0⑴=—1.0062占了)=〃+0618(b—a)=—0.9909

则

(l>

/(孙⑵)=-0.9998>/(x|1)=-1.0000

所以令b=x”⑵=—0.9909,即新的区间为[-1.0062,-0.9909]

区间长度为：1-0.9909-(-1.0062)1=0.0153<e=0.05

所以，最优解：x*=(-0.9909-1.0062)/2=-0.9986

f(x*)=-0.9986*(2-0.9986)=1

第四章线性规划与二次规划

4-1.用单纯形法求解以下线性规划：

/

FindX=[x1,x2,x3]

Minf(X)=-(x]+x2+x3)

io

(2/3)X2<35

(8/15)x2+(2/5)x3<70

(2/15)x2+(3/5)x3<70

(3/4)x,+(7/12)x2+(1/4)尤3<90

(l/4)x,+(5/12)X2+(3/4)X3<90

%1>0

x2>0

x3>0

解：引入非负松弛变量乙,化不等式约束为等式约束:

2

yx2+x4=35

82”

X2+->^3+=70

23

%2+"^"工3+%6=70

371nA

-XH-------H—心+X—90

4},1224377

1532

一Xi4-------Xj~\-----XR4~Xo—90

41122438

构造单纯形表:

X?di

X,X3X4X5X6X7X8

(1)02/301000035

(2)08/152/50100070

(3)02/153/50010070

(4)3/47/121/40001090

(5)1/45/123/40000190

(6)-1-1-1000000

第一轮迭代：

（1）、检查目标行，1、2、3列中负数均为-1,选第一列进基。

（2）检查第一列，对%>0,计算4/“：乙=120,4/%I=360,第四行对应值最小，因此

（3）、以为轴进行高斯消元得下表:

运算）X7di

XX2X3X4X5X6Xg

(7)(1)02/301000035

(8)(2)08/152/50100070

(9)(3)02/153/50010070

(10)(4)*4/317/91/30004/30120

(H)(5)-(10)*l/402/92/3000-1/3160

11

(12)(6)+(10)0-2/9-2/30004/30120

第二轮迭代：

(I),检查目标行，-2/3最小，p=3,x3进基。

(2)、计算第三列4-3:/。23=175,33/。33=350/3,〃/。43=360,^5/%3=90,由最小者决

定q=5o

(3)、以。53为轴进行高斯消元的新表:

运算

X,X2X3X4X5X6X7Xg

(13)⑺02/301000035

(14)(8)-(17)*2/502/500101/5-3/534

(15)(9)-(17)*3/50-1/1500013/10-9/1016

(16)(10)-(17)*1/312/300003/2-1/290

(17)(11)*3/201/31000-1/23/290

(18)(12)+(17)*27300000011180

目标行已无负非基本变量系数，故已得最优解:

X,*=90,X2*=0,X3*=90,/*=-180

第七章机械结构优化

7-1.为什么说满应力法是感性准则法？

答：1、满应力优化设计方法是对结构布局已定的构件尺寸优化，其目的是使结构体枳最小(重量最轻),

主要是针对杆系结构，尤其是桁架结构。

2、它是从力学概念出发建立的一些准则，而不是从最优性条件出发，没有与优化目标函数直接建

立关系，故不能保证目标函数趋于最优，只是以结构满足了力学准则的近似最优设计，它可能

找的是极值点而非最值点，所以说是感性准则法，不大可靠。

7-2.为什么说导重准则法克服了虚功准则法不可克服的缺陷？

答：1、虚功准则法的特点是结构位移采用虚功表达，认为外载荷不随设计变量变化。这对外载荷包括

自重及惯性载荷的航空航天结构、精密机械结构是不成立的。不能考虑载荷对设计变量的导数，这

是虚功准则法无法克服的缺陷。

虚功准则法有如下局限：

(1)对于航空航天、精密机械等惯性载荷结构，虚功准则法具有先天缺陷，准则不准，最优解

不优。

(2)不能对••般混合结构(含非杆、板单元)进行结构优化。

(3)不能进行几何变量优化，设计变量尤不能是坐标。

(4)不能进行动力特性优化。

2、导重准则法克服了虚功准则法的缺陷失因为：

(1)导重准则法考虑了设计变量变化引起的载荷变化。用于优化对自重等惯性载荷为主的航空航

天结构、精密机械结构，效果显著。可保证求得最优解。

(2)借助于线性互补问题的克莱姆算法求解多个不等式约束的库恩-塔克乘子，自动而有效的区分

12

了临界约束与非临界约束。

(3)设计变量除构件尺寸外，还可包括节点坐标。

导重法则是严密推导的数学准则法，能够满足K—T条件，最优结构按各组构件的导重正比分配结

构重量(导重为引导各组构件重量分配，使结构趋于最优)，尤其是对自重及惯性载荷为主的航空航天

结构，精密机械结构尤其适用。所以说导重准则法克服了虚功准则法的先天不足。导重法可考虑多种性

态函数，位移、应力、安全度、精度、基频等，这是虚功法无法做到的，关键是不用位移的虚功表达，

而采用一般刚度方程表达，放开了手脚，使得导重法具有了普适性、通用性。

7-3.导重的意义是什么？单约束与多约束导重法是如何实现不等式约束的？

答：1、导重的意义：导重起到引导各组构件重量分配，使结构趋于最优化的作用。

2、单约束不等式的约束：是通过控制晨G的正负来实现的。

其中QG为目标函数最优值的灵敏度。最优设计点上的目标函数值即目标函数最优值对约束界的灵敏

度等于该约束的负梯度：

即对约束WWW。约束界是W0，有

GAi=—4名生；总导重G=£GXj

'dAi乙

保重设计时，4G的取值正负不加控制。

优重设计时，对于不等式约束：W(X)-可通过控制4G的正负来实现。

①4G>0时，说明卬(X)=W0增加，这使/(X*)减小，但约束限制卬=%不能

再增加，这时人>0,满足K-T条件。

②4G<0,不满足K-T条件之420,它对应W(X)=W。减少，方可使/(X*)改善的情况。

而原不等式约束W(X)4W0是允许W下降的，即设计点X可离开W=W0约束面，向可行域内的无约

束极值点X**方向移动，直至X*=X**时，4=G=0,这样使K-T条件420得到满足，结构重量小

于给定指标，目标函数反而更优，这可称之为“优重设计

这样Kuhn-Tucker5个条件通过以上迭代控制得到了充分的满足，不等式约束极值条件得以实现。

7-4.导重法的特点、优点是什么？有什么地方应与改进？

答：导重法的特点和优点：

(1)导重准则法克服了虚功准则法的先天不足，尤其是对自重等惯性载荷为主的航空航天结构,

精密机械结构，尤为适用。

(2)导重法是严格推导的数学准则法，可保证求得最优解。导重准则物理意义明确、直观、表

13

达简洁。

（3）在计算工作量上与虚功准则法相比，导重法要计算灵敏度，计算量与设计变量数目微弱有

关，虚功准则法要计算虚内力，计算量与约束数目微弱有关。

（4）导重法可考虑多种性态约束，如位移、应力、安全度、精度、基频等，导重法具有普适性、

通用性。而虚功准则法无法做到。

需要改进的地方：此方法只适用于局部寻优。

第八章结构优化的灵敏度分析

8-1、比较各种情况下敏度载荷法与虚载荷法的计算工作量？

答：1）单位移虚载荷法比全位移敏度载荷法的计算工作量大。因为单位移虚载荷法计算的是单元内的

每个自由度所对应位移的敏度，而全位移敏度载荷法则是计算所有节点位移的敏度。

2）对于刚度敏度、载荷阵的敏度、应力敏度与谐振频率的敏度来说，虚载荷法比敏度载荷法的计

算工作量要小。因为敏度载荷都是在虚载荷的基础上计算而得。

8-2、试用虚载荷法推导全位移对全变量的敏度计算表达式

即叫…?

解：由刚度方程KU=P求导得：

euJdP豕]

---二K---------------U

5Xj（以SXj）

人、dP5K

令：P=------------UTT则:①

5XjdXj

再利用虚载荷法，在①式两端同乘eT=[O...1...0],

其中与Um对应的第m个元素为1,其余元素均为0。可得:

=eT.K-|.—-eT.K-1.—.U,两边转置得:

SXjdx-5Xj

/、i

令―,昨图-噫得:

7=P；'UV...................................................................................②

dx-

令X=[xjx2...XN]NX1PMXN=[P]P2...PN]则:

14

3u..v

=PU.③

ax

因为剪="曳也

axdx、5umAX

把①、②、③带入左式，得:

7

auaudx,dum=K-'P/^-7-P'U=K-'^P

PiUvR

axdx,aumax

HSU1...

=>P.——=K-'P-P

1ax1④

又因为P'MxN=[PjP；P；]，PM*N=B；P；…P」

把上式带入④式，得：

—=K-'P

dX

全位移对全变量的敏度计算表达式为：—=K-'P

ax

15

第二、机械优化设计3道练习题

222

3-2.求函数/(X)=X；-2%1X2-3X1X2+%1+5X2+4.5x1-4x2+5的极小值，初始点为[-2,2]T,

误差£不大于0.001o

解：

1、建立数学模型：

/

findx-[x1,x2]

xx

min/(X)=x：—2x^x2-3再々+\+^2+4.5%一4x2+5

s.t.f<0.001初始点为[-2,2]T

2、优化方法：此问题为无约束非线性优化问题，用鲍威尔法进行迭代求解。

3、目标函数的三维图如下：

4、目标函数等值线图如下:

5、迭代输出结果:

16

山上表可知，经5轮迭代可得最优结果：

XI=-0.896104X2=0.291908Fmin=1.9897

6、迭代曲线图：

从中我们可以看出：通过鲍威尔法迭代5次可得出结果，第五次时，目标函数最优值

Fmin=-1.9897,且满足J40.01。

7、迭代原程序如下：

#include<math.h>

#include<iostream.h>

inti,j.N,M.NUM;

doubleS[100],SS[100][100],X0[100],X[100],Xl[100],X2[100],XK[100];

doubleX3[100],XXI[100],XX2[100],XX3[100],A[100],B[100],XPl[100],XP2[100];

doubleF,F0,Fl,F2,F3,FX,DLT,DF,El,EP,L,SDX,Q,D,H,HO;

17

〃N-维数

〃X0[100]-初始点

〃H0-进退法的初始进退距

〃A[100],B[100]--维搜索中区间端点坐标

//E1-一维迭代精度

〃EP-鲍威尔法的终止迭代精度

voidOBJFX();〃目标函数

voidSEEK();〃进退法一维搜索区间

voidGOLDO;//黄金分割法

intmain(intargc,char*argv口)

(

〃输入初始数据

NUM=-1;

N=2;

cout<〈”请输入初始点："<<endl;

for(i=l;i〈=N;i++)

cin»X0[i];

El=0.001;

EP=0.001;

H0=l;

for(i=l;i<=N;i++)

{

for(j=l;j<=N+l;j++)

SS[i][j]=O;

SS[i][i]=l;

}

begin:

for(i=l;i<=N;i++)

(

XXl[i]=X0[i];

X[i]=XXl[i];

)

NUM=NUM+1;

cout<〈"第〈"轮迭代："<<endl;

cout«，zXl="«X0[l]«，/X2="«X0[2];

OBJFXO;

cout«zzF=z，«FX«endl;

F1=FX;

FO=F1;

DLT=T;

for(j=l;j<=N;j++)

(

for(i=l;i<=N;i++)

S[i]=SS[i][j];

18

XO[i]=X[i]：

)

SEEKO;

GOLD();

OBJFXO;

F=FX;

DF=F0-F;

FO=F;

if(DF>DLT)

(

DLT=DF;

M=j;

)

)

SDX=O;

for(i=l;i<=N;i++)

SDX=SDX+pow((X[i]-XXl[i]),2);

if(sqrt(SDX)>EP)

(

for(i=l;i<=N;i++)

XX2[i]=X[i];

F2=F;

for(i=l;i<=N;i++)

(

SS[i][N+l]=X[i]-XXl[i];

S[i]=SS[i][N+l];

}

SEEKO；

GOLD();

for(i=l;i<=N;i++)

(

XK[i]=X[i];

)

for(i=l;i<=N;i++)

(

XX3[i]=2*XX2[i]-XXl[i];

X[i]=XX3[i];

)

OBJFXO;

F3=FX;

Q=(F1-2*F2+F3)*pow((F1-F2-DLT),2);

D=O.5*DLT*pow((Fl-F3),2);

if(F3<Fl)

19

if(Q>=D)

gotoLOOP;

for(j=M+l;j<=N+l;j++)

(

for(i=l;i<=N;i++)

SS[i][j-l]=SS[i][jh

)

for(i=l;i<=N;i++)

(

XO[i]=XK[iJ;

}

gotobegin;

)

LOOP:

if(F2>=F3)

gotoLOOP1;

for(i=l;i<=N;i++)

XO[i]=XX2[i];

gotobegin;

LOOP1:

for(i=l;i<=N;i++)

X0[i]=XX3[i];

gotobegin;

)

else

{

cout<<"第〈”轮迭代，即最优结果：“<〈endl;

cout«//Xl=/，«X[l]«//X2="<〈X[2];

OBJFXO;

cout«a，Fmin=w<<FX<<endl;

)

return0;

)

voidSEEKO

(

doublefl,f2,f3;

H=H0;

for(i=l;i〈=N;i++)

(

Xl[i]=X0[i];

X[i]=Xl[i];

)

OBJFXO;

fl=FX;

for(i=l;i<=N;i++)

20

X2[i]=Xl[i]+H*S[i];

X[i]=X2[i];

}

OBJFXO;

f2=FX;

if(fl>f2)

gotolabel;

H=-H0;

for(i=l;i〈=N;i++)

(

X3[i]=Xl[i];

}

f3=f1;

label2:

for(i=l;i〈=N;i++)

{

Xl[i]=X2[i];

X2[i>X3[i];

)

fl=f2;

f2=f3;

label:

H=2*H;

for(i=l;i<=N;i++)

{

X3[i]=X2[i]+H*S[i];

X[i]=X3[i];

)

OBJFXO;

f3=FX;

if(f2>=f3)

gotolabel2;

if(H<0)

|

for(i=l;i<=N;i++)

(

A[i]=X3[i];

B[i]=Xl[i];

)

}

else

for(i=l;i<=N;i++)

21

AEi]=Xl[i];

B[i]=X3[i];

}

}

)

voidGOLD()

(

doubleYl,Y2;

for(i=l;i<=N;i++)

(

XPl[i]=A[i]+0.618*(B[i]-A[i]);

X[i]=XPl[i];

}

OBJFXO：

Y1=FX;

for(i=l;i〈=N;i++)

{

XP2[i]=A[i]+0.382*(B[i]-A[i])：

X[i>XP2[i];

)

OBJFXO;

Y2=FX;

do

(

if(YK=Y2)

(

for(i=l;i<=N;i++)

(

A[i]=XP2[i];

XP2[i]=XPl[i];

)

Y2=Y1;

for(i=l;i<=N;i++)

(

XPl[i]=A[i]+0.618*(B[i]-A[i]);

X[i]=XPl[i];

)