最优化-刘志斌-课后习题3-5参考答案
练习题三
1、用0.618法求解问题
min(t)t32t1
t0
的近似最优解,已知(t)的单谷区间为[0,3],要求最后区间精度0.5。 答:t=0.8115;最小值-0.0886.(调用golds.m函数)(见例题讲解5) 2、求无约束非线性规划问题
22
min f(x1,x2,x3)=x124x2x32x1
的最优解
解一:由极值存在的必要条件求出稳定点:
fff
2x12,8x2,2x3,则由fx0得x11,x20,x30 x1x2x3
再用充分条件进行检验:
2f2f2f2f2f2f
0,0 2,28,22,0,
x1x3x2x3x12x2x3x1x2
200
即2f080为正定矩阵得极小点为x*(1,0,0)T,最优值为-1。
002
解二:目标函数改写成
22
min f(x1,x2,x3)=(x11)24x2x31
易知最优解为(1,0,0),最优值为-1。
3、用最速下降法求解无约束非线性规划问题。
2
minf(X)x1x22x122x1x2x2
其中X(x1,x2)T,给定初始点X0(0,0)T。
f(x)
(x)14x2x
112
解一:目标函数f(x)的梯度f(x) f(x)12x12x2(x)2
11
f(X(0))令搜索方向d(1)f(X(0))再从X(0)出发,沿d(1)方向作一维寻
1101(0)(1)
优,令步长变量为,最优步长为1,则有Xd
01故f(x)f(X(0)d(1))()2()22()2221()
011
令1'()220可得11 X(1)X(0)1d(1) 求出X(1)点之后,
01111
与上类似地,进行第二次迭代:f(X(1)) 令d(2)f(X(1))
11令步长变量为,最优步长为2,则有
X故
(1)
d
(2)
111
111
f(x)f(X(1)d(2))(1)(1)2(1)22(1)(1)(1)252212()0.8111(2)(1)(2)1令()1020可得 2 XX2d 5111.25
'
2
0.2(2)
f(X)0.2828 此时所达到的精度f(X(2))
0.21
本题最优解X,f(X)1,25
1.5
解二:利用matlab程序求解
首先建立目标函数及其梯度函数的M文件 function f=fun(x)
f=x(1)-x(2)+2*x(1)*x(1)+2*x(1)*x(2)+x(2)*x(2); function g=gfun(x)
g=[1+4*x(1)+2*x(2),-1+2*x(1) +2* x(2) ]; 调用grad.m文件 x0=[0,0];
[x,val,k]=grad('fun','gfun',x0)
结果
x=[ -1.0000 ,1.5000] val= -1.2500 k=33
即迭代33次的到最优解x=[ -1.0000 ,1.5000];最优值val= -1.2500。
4、试用Newton法求解第3题。 解一:计算目标函数的梯度和Hesse阵
f(x)
(x)14x2x
112
目标函数f(x)的梯度f(x) f(x)12x12x2(x)2420.50.51
,其逆矩阵为 f(X)GG1220.5
2
0.50.5TTT
X(1)X(0)G1f(X(0))0,01,11,1.50.51计算f(X(1))0。
1
本题最优解X,f(X)1,25
1.5
解二:除了第3题建立两个M文件外,还需建立Hesse矩阵的M文件 利用matlab程序求解
首先建立目标函数及其梯度函数的M文件 function f=fun(x)
f=x(1)-x(2)+2*x(1)*x(1)+2*x(1)*x(2)+x(2)*x(2); function g=gfun(x)
g=[1+4*x(1)+2*x(2),-1+2*x(1) +2* x(2) ]; function h=hess(x) g=[4 2;2 2 ]; 调用newton.m文件
x0=[0,0];
[x,val,k]=newton('fun','gfun','hess',x0) 结果
x=[ -1.0000 ,1.5000] val= -1.2500 k=1
5、用Fletcher—Reeves法求解问题
2 minf(X)x1225x2
其中X(x1,x2)T,要求选取初始点X0(2,2)T,106。 解一:
20x1201
,rf(x)(2x1,50x2)T. f(x)(x1,x2),G2050050x2第一次迭代:令p0r0(4,100)T,
4
(4,100)100
r0Tr010T p0Gp020450
(4,100)
050100即,X(1)X(0)0p0(1.92,0)T 第二次迭代:
||r1||23
,p1r10p0(3.846,0.15)T r1(3.84,0),02
||r0||2000
T
3.84(3.84,0)
0r1Tr11T0.4802
p1Gp1203.846
(3.846,0.15)
0500.15
X(2)X(1)1p1(0.0732,0.072)T
第三次迭代:
r2(0.1464,3.6)T……(建议同学们自己做下去,注意判别rk)
解二:利用matlab程序求解
首先建立目标函数及其梯度函数的M文件 function f=fun(x) f=x(1)^2+25* x(2)*x(2); function g=gfun(x) g=[2*x(1), 50* x(2) ]; 调用frcg.m文件 x0=[2,2]’;epsilon=1e-6;
[x,val,k]=frcg('fun','gfun',x0, epsilon) 结果
x = 1.0e-006 *[ 0.2651, 0.0088] val =7.2182e-014 k = 61
6、试用外点法(二次罚函数方法)求解非线性规划问题
2
minf(X)(x12)2x2 s.t.g(X)x210
其中X(x1,x2)R2
解:设计罚函数 P(x,M)
f(X)
M*[g(X) ]
采用Matlab编程计算,结果x=[1 0];最优结果为1。(调用waidianfa.m) 7、用内点法(内点障碍罚函数法)求解非线性规划问题:
min(x11)3x2
tx110s.. x20
解:容易看出此问题最优解为x=[1 0];最优值为8.
)1(x31)2xr(1x(给出罚函数为 P(x,r1/
1) 2x1/
)
令
P(x,r)rP(x,r)r
3(x11)2010 ;2
x1x2(x11)2x2
1/2
1(1
从而当r
0时,x(r)x
0
(建议同学自己编程序计算)
22
minf(X)x1x2
8、用乘子法求解下列问题
h(X)xx20112
解:建立乘子法的增广目标函数:
2
(x,,)x12x2(x1x22)
2
(x1x22)^2
令:
(x,,)
2x1(x1x22)0 x1
(x,,)
2x2(x1x22)0 x1
解上述关于x的二元一次方程组得到稳定点
222
222
1
当乘子取2时,或发参数趋于无穷时,得到x*即最优解。
1(建议同学自己编程序计算)
练习题四
1、石油输送管道铺设最优方案的选择问题:考察网络图4-6,设A为出发地,F为目的地,B,C,D,E分别为四个必须建立油泵加压站的地区。图中的线段表示管道可铺设的位置,线段旁的数字表示铺设这些管线所需的费用。问如何铺设管道才能使总费用最小?
图4- 1
解:
第五阶段:E1—F 4;E2—F 3;
第四阶段:D1—E1 —F 7;D2—E2—F 5;D3—E1—F 5;
第三阶段:C1—D1—E1 —F 12;C2—D2—E2—F 10;C3—D2—E2—F 8;C4—
D3—E1—F 9;
第二阶段:B1—C2—D2—E2—F 13; B2—C3—D2—E2—F 15; 第一阶段:A—B1—C2—D2—E2—F 17; 最优解:A—B1—C2—D2—E2—F 最优值:17
2、 用动态规划方法求解非线性规划
maxf(x)x1x2x327
x1,x2,x30
解:x19,x29,x39,最优值为9。 3、用动态规划方法求解非线性规划
2
maxz7x126x15x2
tx12x210s..
x3x912
x1,x20
解:用顺序算法
阶段:分成两个阶段,且阶段1 、2 分别对应x1,x2。 决策变量:x1,x2
状态变量:vi,wi分别为第j 阶段第一、第二约束条件可供分配的右段数值。
f1*(v1,w1)max{7x126x1}min{7v126v1,7w126w1}
0x1v1
0x1w1
*x1min{v1,w1}
2f2*(v2,w2)max{5x2f1*(v22x2,w23x2)}
0x25
max{5xmin{7(v22x2)6(v22x2),7(w23x2)6(w23x2)}}
0x25
2
2
22
由于v210,w29,
22
f2*(v2,w2)f2*(10,9)max{min{33x2292x2760,68x2396x2621}
0x25
可解的x19.6,x20.2,最优值为702.92。
4、设四个城市之间的公路网如图4-7。两点连线旁的数字表示两地间的距离。使用迭代法求各地到城市4的最短路线及相应的最短距离。
27
58
61
3
4
图4- 2 城市公路网
解:城市1到城市4路线——1-3-4 距离10;
城市2到城市4路线——2-4 距离8;城市3到城市4路线——3-4 距离4。 5、某公司打算在3个不同的地区设置4个销售点,根据市场部门估计,在不同地区设置不同数量的销售点每月可得到的利润如表4-19所示。试问在各地区如何设置销售点可使每月总利润最大。
表4- 1
解:
将问题分为3个阶段,k=1,2,3;
决策变量xk表示分配给第k个地区的销售点数;
状态变量为sk表示分配给第k个至第3个地区的销售点总数; 状态转移方程:sk+1=sk-xk,其中s1=4; 允许决策集合:Dk(sk)={xk|0≤xk≤sk}
阶段指标函数:gk(xk)表示xk个销售点分配给第k个地区所获得的利润;
最优指标函数fk(sk)表示将数量为sk的销售点分配给第k个至第3个地区所得到的最大利润,动态规划基本方程为:
[gk(xk)fk1(skxk)] k3,2,1fk(sk)0maxxksk
f4(s4)0
k=3时,f3(s3)max[g3(x3)]
x3s3
g3(x3)
00
0x2s2
110
234
f3(s3) x3*
14
16
17
001
234
k=2时,f2(s2)max[g2(x2)1f3(s2x2)]
2
34
10141617
g2(x2)+f3(s2-x2)
01234
000
k=1时,f1(s1)max[g1(x1)f2(s1x1)],f1(s1)max[g1(x1)f2(4x1)]
s0x4
10x0+1012+012
1
1
1
f2(s2) x2*
112
234
0+1412+1017+00+1612+1417+1021+00+17 12+16 17+14 21+10 22+0
222731
2,3
最优解为:x1*=2,x2*=1,x3*=1,f1(4)=47,即在第1个地区设置2个销售点,第2个地区设置1个销售点,第3个地区设置1个销售点,每月可获利润47。
6、设某厂计划全年生产某种产品A。其四个季度的订货量分别为600公斤,700
公斤,500公斤和1200公斤。已知生产产品A的生产费用与产品的平方成正比,系数为0.005。厂内有仓库可存放产品,存储费为每公斤每季度1元。求最佳的生产安排使年总成本最小。
解:四个季度为四个阶段,采用阶段编号与季度顺序一致。 设 sk 为第k季初的库存量,则边界条件为 s1=s5=0
设 xk 为第k季的生产量,设 yk 为第k季的订货量;sk ,xk ,yk 都取实数,状态转移方程为 sk+1=sk+xk - yk 仍采用反向递推,但注意阶段编号是正向的目标函数为:
f1(x)min
x1,x2,x3,x4
(0.005x
i1
4
2
i
si)
第一步:(第四季度) 总效果 f4(s4,x4)=0.005 x42+s4
由边界条件有: s5= s4 + x4 – y4=0,解得:x4*=1200 – s4 将x4*代入 f4(s4,x4)得:
f4*(s4)=0.005(1200 – s4)2+s4=7200 –11 s4+0.005 s42 第二步:(第三、四季度) 总效果 f3(s3,x3)=0.005 x32+s3+ f4*(s4) 将 s4= s3 + x3 – 500 代入 f3(s3,x3) 得:
2
f3(s3,x3)0.005x3s3720011(x3s3500)
0.005(x3s3500)2
22
0.01x30.01x3s316x30.005s315s313950
f3(s3,x3)
0.02x30.01s3160x3解得
x38000.5s3,
代入f3(s3,x3)得
2
f3(s3)75507s30.0025s3
第三步:(第二、三、四季度) 总效果 f2(s2,x2)=0.005 x22+s2+ f3*(s3) 将 s3= s2 + x2 700 代入 f2(s2,x2) 得:
2
f2(s2,x2)0.005x2s275507(x2s2700)
0.0025(x2s2700)2
f2(s2,x2)
0.015x20.005(s2700)70 x2解得
x2700s2,
代入f2(s2,x2)得
2
f2(s2)100006s23)s2
第四步:(第一、二、三、四季度) 总效果 f1(s1,x1)=0.005 x12+s1+ f2*(s2)
将 s2= s1 + x1 – 600= x1 – 600 代入 f1(s1,x1) 得:
f1(s1,x1)0.005x12s1100006(x1600)
3)(x1600)2
f1(s1,x1)
(0.043)x180x1解得
x1600,
代入f1(s1,x1)得
f1(s2)11800
由此回溯:得最优生产–库存方案
x1*=600,s2*=0; x2*=700,s3*=0; x3*=800,s4*=300; x4*=900。
7、某种机器可在高低两种不同的负荷下进行生产。设机器在高负荷下生产的产量函数为g=8u1,其中u1为投入生产的机器数量,年完好率a=0.7;在低负荷下生产的产量函数为h=5y,其中y为投入生产的机器数量,年完好率为b=0.9。假定开始生产时完好机器的数量s1=1000。试问每年如何安排机器在高、低负荷下的生产,使在5年内生产的产品总产量最高。 解:
构造这个问题的动态规划模型: 设阶段序数k表示年度。
状态变量sk为第k年度初拥有的完好机器数量,同时也是第k−1年度末时的完好机器数量。
决策变量uk为第k年度中分配高负荷下生产的机器数量,于是sk−uk为该年度中分配
在低负荷下生产的机器数量。
这里sk和uk均取连续变量,它们的非整数值可以这样理解,如sk=0.6,就表示一台机器在k年度中正常工作时间只占6/10;uk=0.3,就表示一台机器在该年度只有3/10的时间能在高负荷下工作。
状态转移方程为:sk1aukb(skuk)0.7uk0.9(skuk), k1,2,,5 k段允许决策集合为:Dk(sk)uk0uksk 设vk(sk,uk)为第k年度的产量,则vk8uk5(skuk) 故指标函数为:V1,5vk(sk,uk)
k1
5
令最优值函数fk(sk)表示由资源量sk出发,从第k年开始到第5年结束时所生产的产品的总产量最大值。因而有逆推关系式:
f6(s6)0
8uk5(skuk)fk10.7uk0.9(skuk) fk(sk)umaxkDk(sk)
k1,2,3,4,5
从第5年度开始,向前逆推计算。 当k=5时,有:
f5(s5)max8u55(s5u5)f60.7u50.9(s5u5)
0u5s50u5s50u5s5
max8u55(s5u55)max3u55s5
因f5是u5的线性单调增函数,故得最大解u5*,相应的有:
f5(s5)8s5
当k=4时,有:
f4(s4)max8u45(s4u4)f50.7u40.9(s4u4)
0u4s40u4s40u4s40u4s4
max8u45(s4u4)80.7u40.9(s4u4)max13.6u412.2(s4u4)max1.4u412.2s4
故得最大解,u4*=s4,相应的有
f4(s4)13.6s4
依此类推,可求得
*u3s3, 相应的 f3(s3)17.5s3*
u20, 相应的 f2(s2)20.8s2 *
u10, 相应的 f1(s1)23.7s1
因s1=1000,故:f1(s1)23700 计算结果表明:最优策略为
*****u10,u20,u3s3,u4s4,u5s5
即前两年应把年初全部完好机器投入低负荷生产,后三年应把年初全部完好机器投入高负荷生产。这样所得的产量最高,其最高产量为23700台。
在得到整个问题的最优指标函数值和最优策略后,还需反过来确定每年年初的状态,即从始端向终端递推计算出每年年初完好机器数。已知s1=1000台,于是可得:
**
s20.7u10.9(s1u1)0.9s1900(台)**s30.7u20.2(s2u2)0.9s2810(台)**s40.7u30.9(s3u3)0.7s3567(台) **s50.7u40.9(s4u4)0.7s4397(台)**s60.7u50.9(s5u5)0.7s5278(台)
8、有一辆最大货运量为10t 的卡车,用以装载3种货物,每种货物的单位重量及相应单位价值如表4-20所示。应如何装载可使总价值最大?
表4- 2
解:建模设三种物品各装x1,x2,x3件
max(4x15x26x3)3x14x25x310
xj0,xjI,j1,2,3
利用动态规划的逆序解法求此问题。
s1c,D1(s1){x1|0x1s1} s2s1x1,D2(s2){x2|0x2s2} s3s2x2,D3(s3){x3|0x3s3}
状态转移方程为: sk1Tk(sk,xk)skx,kk3, 2,1
该题是三阶段决策过程,故可假想存在第四个阶段,而x40,于是动态规划的基本方程为:
[xk,fk1(sk1)],k3,2,1fk(sk)xkmaxDK(sk)
f4(s4)0
k3,
f3(s3)k2,
x30,1,,[s3/5]
max
6x3
f2(s2)
k1,
max[5x2f3(s3)]
x20,1,,[
s2
]4
max[5x2f3(s24x2)]
x20,1,,[
s2]4
f1(s1)max[4x1f2(s2)]
x10,1,2,3
max[4x1f2(s13x1)]
x10,1,2,3
计算最终结果为x12,x21,x30,最大价值为13 。
9、设有 A,B,C 三部机器串联生产某种产品,由于工艺技术问题,产品常出现次品。统计结果表明,机器 A,B,C产生次品的概率分别为 pA=30%, PB=40%, PC=20%, 而产品必须经过三部机器顺序加工才能完成。为了降低产品的次品率,决定拨款 5 万元进行技术改造,以便最大限度地提高产品的成品率指标。现提出如下四种改进方案:
方案1:不拨款,机器保持原状;
方案2:加装监视设备,每部机器需款 1 万元; 方案3:加装设备,每部机器需款 2 万元;
方案4:同时加装监视及控制设备,每部机器需款 3 万元;
采用各方案后,各部机器的次品率如表4-21。
表4- 3
问如何配置拨款才能使串联系统的可靠性最大?
解:为三台机器分配改造拨款,设拨款顺序为A, B, C,阶段序号反向编号为 k,即第一阶段计算给机器 C 拨款的效果。
设 sk 为第 k 阶段剩余款,则边界条件为 s3=5; 设 xk 为第 k 阶段的拨款额; 状态转移方程为 sk-1=sk-xk;
目标函数为 max R=(1-PA)(1-PB)(1-PC) 仍采用反向递推 第一阶段 :对机器 C 拨款的效果
R1(s1,x1)=d1(s1,x1) R0(s0,x0)= d1(s1,x1)
第二阶段 :对机器 B, C 拨款的效果 由于机器 A 最多只需 3 万元,故 s2 2 递推公式:
R2(s2,x2)=d2(s2,x2) R1(s1,x1*)
例:R2(3,2)=d2(3,2) R1(1,1)=(1-0.2) 0.9=0.72 得第二阶段最优决策表
第三阶段 :对机器 A, B, C 拨款的效果 边界条件:s3 = 5 递推公式:
R3(s3,x3)=d3(s3,x3) R2(s2,x2*)
例:R3(5,3)=d3(5,3) R2(2,2)=(1-0.05) 0.64=0.608 得第三阶段最优决策表
回溯 :有多组最优解。
I:x3=1, x2=3, x1=1, R3=0.8 0.9 0.9=0.648 II:x3=2, x2=2, x1=1, R3= 0.90.80.9=0.648 III: x3=2, x2=3, x1=0, R3= 0.90.90.8=0.648
练习题五
1、考察多目标规划问题
x2,2x1
1,1x2,试画出个目标函数的图形,并求出其中f1(x)x2,f2(x)
x1,x2
R1,R2,Rab,Rpa,Rwp,这里Ri是minfi(x)的最优解集。
2x4
解:
2、用线性加权法中的法求解下述多目标规划问题
minf1(x)4x16x2maxf2(x)3x13x2。 2x14x214s..t6x13x224x,x012
解:minf1(x)4x16x2最优解为x(1)=[0 0]T;
maxf2(x)3x13x2最优解为x(2)=[3 2]T; 利用法得线性方程组:
1*02*0
1*242*15 112
解得唯一加权系数[0.3846,0.6154] 原多目标规划加权后
T
minF(x)0.3846f1(x)0.6154f2(x)
解得加权后的最优解为:x=[4 0]T,最优值为-1.2312
3、用线性加权求和法求解下述多目标规划问题,取10.6,20.4。
vmin
F(x)(x13x2,2x1x2)T
3x12x26x3x32s..t1
2x1x22x1,x20
解:将问题转化为一个新的单目标规划问题。
minv(F(x))0.6(x13x2)0.4(2x1x2)
约束条件同上,该问题转化为线性规划问题,可用单纯形法求解,也可用Matlab命令求解(求解过程略)。
解得加权后的最优解为:x=[0 1]T,最优值为-1.4。 4、用平方和加权法求解多目标规划问题:
T
Vmin(f(x),f(x))12xD
x1x24
12
其中 f1(x)x1,f2(x)x2,D:x1x28,1,2。
33x,x0
12
解:不难看出两个目标函数下界均为0,得平方和加权法后的新目标规划问题:
122
minF(x)x12x2
33
x1x24
D:x1x28 x,x012
利用matlab程序求解
首先建立目标函数及其梯度函数的M文件 function f=fun(x)
f=1/3*x(1)^2+2/3* x(2)*x(2);
[x,fval]=fmincon(‘f’,[0 0],[1 -1;1 1],[4;8],[],[],[0 0]) 解得最优解为:x=[0 0]T,最优值为0。
5、用极小极大法和Matlab软件求解下述多目标规划问题
vmins..t
2
F(x)((x13)2x2,x12(x22)2)T
x1x22
i
。
2
,x12(x22)2],再求 解:取评价函数为v(F(x))max[(x13)2x2
2
minv(F(x))min{max[(x13)2x2,x12(x22)2]}
D
i
Matlab软件求解: 编制M文件 function f=mnmax(x) f(1)=(x(1)-3)^2+x(2)^2; f(2)=x(1)^2+(x(2)-2)^2 设初值 x0=[0;0]; 调用函数
[x,fval]=fminimax(@mnmax,x0,[1 1],[2]) 结果: x = 1.3000
0.7000 fval =
3.3800 3.3800
可得X[1.3,0.7]T;对应f1f23.38从而X[1.3,0.7]T为原问题的解。
附习题中用过的Matlab程序
1、bbmethod
function [x,y]=bbmethod(f,G,h,Geq,heq,lb,ub,x,id,options)
%整数线性规划分支定界法,可求解纯整数规划和混合整数规划。
%y=minf’*x s.t. G*x
%[x,y]=bbmethod(f,G,h,Geq,heq,lb,ub,x,id,options)
%参数说明
%lb:解的下界列向量(Default:-int)
%ub:解的上界列向量(Default:int)
%x:迭代初值列向量
%id:整数变量指标列向量,1-整数,0-实数(Default:1)
global upper opt c x0 A b Aeq beq ID options;
if nargin
options.LargeScale='off';end
if nargin
if nargin
if nargin
if nargin
if nargin
if nargin
upper=inf;c=f;A=G; b=h;Aeq=Geq;beq=heq;x0=x;ID=id;
ftemp=IntLP(lb(:),ub(:));
x=opt;y=upper;
%下面是子函数
function ftemp=IntLP(vlb,vub)
global upper opt c x0 A b Aeq beq ID options;
[x,ftemp,how]=linprog(c,A,b,Aeq,beq,vlb,vub,x0,options);
if how
return;
end;
if ftemp-upper>0.00005 %in order to avoid error
return;
end;
if max(abs(x.*ID-round(x.*ID)))
if upper-ftemp>0.00005 %in order to avoid error
opt=x';upper=ftemp;
return;
else
opt=[opt;x'];
return;
end;
end;
notintx=find(abs(x-round(x))>=0.00005); %in order to avoid error
intx=fix(x);tempvlb=vlb;tempvub=vub;
if vub(notintx(1,1),1)>=intx(notintx(1,1),1)+1;
tempvlb(notintx(1,1),1)=intx(notintx(1,1),1)+1;
ftemp=IntLP(tempvlb,vub);
end;
if vlb(notintx(1,1),1)
tempvub(notintx(1,1),1)=intx(notintx(1,1),1);
ftemp=IntLP(vlb,tempvub);
end;
2、golds.m
function [s,phis,k,G,E]=golds(phi,a,b,delta,epsilon)
%功能: 0.618法精确线搜索
%输入: phi是目标函数, a, b 是搜索区间的两个端点
% delta, epsilon分别是自变量和函数值的容许误差
%输出: s, phis分别是近似极小点和极小值, G是nx4矩阵,
% 其第k行分别是a,p,q,b的第k次迭代值[ak,pk,qk,bk],
% E=[ds,dphi], 分别是s和phis的误差限.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=(sqrt(5)-1)/2; h=b-a; phia=feval(phi,a); phib=feval(phi,b);
p=a+(1-t)*h; q=a+t*h; phip=feval(phi,p); phiq=feval(phi,q);
k=1; G(k,:)=[a, p, q, b];
while(abs(phib-phia)>epsilon)|(h>delta)
if(phip
b=q; phib=phiq; q=p; phiq=phip;
h=b-a; p=a+(1-t)*h; phip=feval(phi,p);
else
a=p; phia=phip; p=q; phip=phiq;
h=b-a; q=a+t*h; phiq=feval(phi,q);
end
k=k+1; G(k,:)=[a, p, q, b];
end
ds=abs(b-a); dphi=abs(phib-phia);
if(phip
s=p; phis=phip;
else
s=q; phis=phiq;
end
E=[ds,dphi];
3、grad.m
function [x,val,k]=grad(fun,gfun,x0)
% 功能: 用最速下降法求解无约束问题: min f(x)
%输入: x0是初始点, fun, gfun分别是目标函数和梯度
%输出: x, val分别是近似最优点和最优值, k是迭代次数.
maxk=5000; %最大迭代次数
rho=0.5;sigma=0.4;
k=0; epsilon=1e-5;
while(k
g=feval(gfun,x0); %计算梯度
d=-g; %计算搜索方向
if(norm(d)
m=0; mk=0;
while(m
if(feval(fun,x0+rho^m*d)
mk=m; break;
end
m=m+1;
end
x0=x0+rho^mk*d;
k=k+1;
end
x=x0;
val=feval(fun,x0);
4、newton.m
function [x,val,k]=newton(fun,gfun,Hess,x0)
%功能: 用尼牛顿法求解无约束问题: min f(x)
%输入: x0是初始点, fun, gfun, Hess 分别是求
% 目标函数,梯度,Hesse 阵的函数
%输出: x, val分别是近似最优点和最优值, k是迭代次数.
maxk=100; %给出最大迭代次数
sigma=0.4;
k=0; epsilon=1e-5;
while(k
gk=feval(gfun,x0); %计算梯度
Gk=feval(Hess,x0); %计算Hesse阵
dk=-Gk\gk'; %解方程组Gk*dk=-gk, 计算搜索方向
if(norm(gk)
x0=x0+dk';
k=k+1;
end
x=x0;
val=feval(fun,x);
5、frcg.m
function [x,val,k]=frcg(fun,gfun,x0)
% 功能: 用FR共轭梯度法求解无约束问题: min f(x)
%输入: x0是初始点, fun, gfun分别是目标函数和梯度
%输出: x, val分别是近似最优点和最优值, k是迭代次数.
maxk=5000; %最大迭代次数
rho=0.6;sigma=0.4;
k=0; epsilon=1e-6;
n=length(x0);
while(k
g=feval(gfun,x0); %计算梯度
itern=k-(n+1)*floor(k/(n+1));
itern=itern+1;
%计算搜索方向
if(itern==1)
d=-g;
else
beta=(g'*g)/(g0'*g0);
d=-g+beta*d0; gd=g'*d;
if(gd>=0.0)
d=-g;
end
end
if(norm(g)
m=0; mk=0;
while(m
if(feval(fun,x0+rho^m*d)
mk=m; break;
end
m=m+1;
end
x0=x0+rho^mk*d;
val=feval(fun,x0);
g0=g; d0=d;
k=k+1;
end
x=x0;
val=feval(fun,x);
6、waidianfa.m
close all
clear all
clc
syms x1 x2 M; %M为罚因子。
m(1)=1;
c=8; %c为递增系数。赋初值。
a(1)=20;
b(1)=20;
f=(x1-2)^2+x2^2+M*(x1-1)^2; %外点罚函数
f0(1)=500;
%求偏导、Hessian元素
fx1=diff(f,'x1');
fx2=diff(f,'x2');
fx1x1=diff(fx1,'x1');
fx1x2=diff(fx1,'x2');
fx2x1=diff(fx2,'x1');
fx2x2=diff(fx2,'x2');
%外点法M迭代循环
for k=1:100 x1=a(k);
x2=b(k);
M=m(k); %牛顿法求最优值
for n=1:100 f1=subs(fx1); %求解梯度值和Hessian矩阵
f2=subs(fx2);
f11=subs(fx1x1);
f12=subs(fx1x2);
f21=subs(fx2x1);
f22=subs(fx2x2);
if(double(sqrt(f1^2+f2^2))
a(k+1)=double(x1);b(k+1)=double(x2);f0(k+1)=double(subs(f));
break;
else
X=[x1 x2]'-inv([f11 f12;f21 f22])*[f1 f2]';
x1=X(1,1);x2=X(2,1);
end
end
if(double(sqrt((a(k+1)-a(k))^2+(b(k+1)-b(k))^2))
%输出最优点坐标,罚因子迭代次数,最优值
X=[a(k+1) b(k+1)]
step=k
F=f0(k+1)
break;
else
m(k+1)=c*m(k);
end
end