共轭梯度法编程.docx
《共轭梯度法编程.docx》由会员分享,可在线阅读,更多相关《共轭梯度法编程.docx(20页珍藏版)》请在冰豆网上搜索。
![共轭梯度法编程.docx](https://file1.bdocx.com/fileroot1/2023-2/6/c6a5bc16-f4da-4a18-9f79-4dd2edd99a45/c6a5bc16-f4da-4a18-9f79-4dd2edd99a451.gif)
共轭梯度法编程
P90第二章11
(2)
用大M法求解
minw=2x1+x2-x3-x4
s.tx1-x2+2x3-x4=2
2x1+x2-3x3+x4=6
x1+x2+x3+x4=7
xi≥0,i=1,2,3,4
用matlab求解如下:
f=[2,1,-1,-1,200,200,200];
a=[1-12-1100;21-31010;1111001];
b=[267];
lb=zeros(7,1);
[x,fval,exitflag,output,lambda]=linprog(f,[],[],a,b,lb)
Optimizationterminated.
运行结果如下:
x=
3.0000
0.0000
1.0000
3.0000
0.0000
0.0000
0.0000
fval=
2.0000
exitflag=
1
output=
iterations:
7
algorithm:
'large-scale:
interiorpoint'
cgiterations:
0
message:
'Optimizationterminated.'
lambda=
ineqlin:
[0x1double]
eqlin:
[3x1double]
upper:
[7x1double]
lower:
[7x1double]
从上述运行结果可以得出:
最优解为x=,最小值约为f*=2。
P151第三章26
用共轭梯度算法求f(x)=(x1-1)^2+5*(x2-x1^2)^2的极小点,取初始点x0=。
用matlab求解如下:
functionmg=MG()
%
%共轭梯度法求解习题三第26题
%
clc;
clear;
n=2;
x=[20]';
max_k=100;
count_k=1;
trace(1,1)=x
(1);
trace(2,1)=x
(2);
trace(3,1)=f_fun(x);
k=0;
g1=f_dfun(x);
s=-g1;
whilecount_k<=max_k
ifk==n
g0=f_dfun(x);
s=-g0;
k=0;
else
r_min=fminbnd(@(t)f_fun(x+t*s),-100,100);
x=x+r_min*s;
g0=g1;
g1=f_dfun(x);
ifnorm(g1)<10^(-6)
break;
end
m=(norm(g1)^2)/(norm(g0)^2);
s=-g1+m*s;
count_k=count_k+1;
trace(1,count_k)=x
(1);
trace(2,count_k)=x
(2);
trace(3,count_k)=f_fun(x);
k=k+1;
end
end
count_k
x
f=f_fun(x)
functiong=f_dfun(x)
g(1,1)=20*x
(1)^3-20*x
(1)*x
(2)+2*x
(1)-2;
g(2,1)=10*x
(2)-10*x
(1)^2;
functionf=f_fun(x)
f=5*x
(1)^4-10*x
(1)^2*x
(2)+x
(1)^2-2*x
(1)+1+5*x
(2)^2;
运行结果如下:
k=
59
x0=
1.0000
1.0000
f0=
0
从上述运行结果可以得出:
最优解为x=,最小值约为f*=0。
P151第三章26
用BFGS算法求f(x)=(x1-1)^2+5*(x2-x1^2)^2的极小点,取初始点x0=。
用matlab求解如下:
symsx1x2;
f=(x1-1)^2+5*(x2-x1^2)^2;
v=[x1,x2];
df=jacobian(f,v);
df=df.';
x0=[2,0]';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});k=0;
H0=[1,0;0,1];
while(norm(g1)>0.001&k<300)
ifk==0
p=-H0*g1;
else
vk=sk/(sk'*yk)-(H0*yk)/(yk'*H0*yk);
w1=(yk'*H0*yk)*vk*vk';
H1=H0-(H0*yk*yk'*H0)/(yk'*H0*yk)+(sk*sk')/(sk'*yk)+w1;
p=-H1*g1;
H0=H1;
end
x00=x0;
result=Usearch1(f,x1,x2,df,x0,p);
arf=result
(1);
x0=x0+arf*p;
g0=g1;
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
p0=p;
yk=g1-g0;sk=x0-x00;
k=k+1;
end;
k
x0
f0=subs(f,{x1,x2},{x0(1,1),x0(2,1)})
functionresult=Usearch1(f,x1,x2,df,x0,p)
mu=0.001;sgma=0.99;a=0;b=inf;arf=1;
pk=p;
x3=x0;
x4=x3+arf*pk;
f1=subs(f,{x1,x2},{x3(1,1),x3(2,1)});
f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});
gk1=subs(df,{x1,x2},{x3(1,1),x3(2,1)});
gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});
while(f1-f2<=-mu*arf*gk1'*pk)
b=arf;arf=(a+arf)/2;x4=x3+arf*pk;f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});
end;
while(1>0)
if(gk2'*pka=arf;a=min(2*arf,(a+b)/2);x4=x3+arf*pk;f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});
while(f1-f2<=-mu*arf*gk1'*pk)
b=arf;arf=(a+arf)/2;x4=x3+arf*pk;f2=subs(f,{x1,x2},{x4(1,1),x4(2,1)});gk2=subs(df,{x1,x2},{x4(1,1),x4(2,1)});
end;
else
break;
end
end;
result=[arf];
运行结果如下:
k=
19
x0=
1.0000
1.0000
f0=
4.4505e-013
从上述运行结果可以得出:
最优解为x=,最小值约为f*=0。
P229第四章10
用SUMT内点法求解
(1)minf(x)=x1+x2
s.tg1(x)=-x1^2+x2≥0
g2(x)=x1≥0
用matlab求解如下:
functionsumt=SUMT(x0,e0,max_k0)
%SUMT内点法
globalxsM
x=x0;
e=e0;
max_k=max_k0;
trace(1,1)=x
(1);
trace(2,1)=x
(2);
M=100;
c=3;
e_FR=10^-10;
max_FR=200;
fork=0:
max_k
x=FR(x,e_FR,max_FR);
trace(1,k+2)=x
(1);
trace(2,k+2)=x
(2);
iff_pfun(x)break;
end
M=c*M;
end
x
f=f_fun(x)
k
trace
functionf=FR(x0,e,max_k)
globalxs;
count_k=1;
k=0;
x=x0;
g1=f_dfun(x);
s=-g1;
whilecount_k<=max_k
ifk==n
g0=f_dfun(x);
s=-g0;
k=0;
else
r_min=fminbnd(@(t)f_fun(x+t*s),-100,100);
x=x+r_min*s;%
g0=g1;
g1=f_dfun(x)
ifnorm(g1)break;
end
m=(norm(g1)^2)/(norm(g0)^2);
s=-g1+m*s;
count_k=count_k+1;
k=k+1;
end
end
count_k=count_k-1
f=x
注意:
以下三个函数体在运行计算习题四第8
(1)题时使用
%*****************************************************
functionf=f_fun(x)
globalM;
f=x
(1)+x
(2)+M*f_pfun(x);%题8
(1)
functiong=f_dfun(x)
globalM;
a=-x
(1)^2+x
(2);
b=x
(1);
ifa>=0&b>=0
g(1,1)=1;
g(2,1)=1;
elseifa<0&b>=0
g(1,1)=1+2*M*(-x
(1)^2+x
(2))*(-2*x
(1));
g(2,1)=1+2*M*(-x
(1)^2+x
(2));
elseifa>=0&b<0
g(1,1)=1+2*M*x
(1);
g(2,1)=1;
else
g(1,1)=1+2*M*(-x
(1)^2+x
(2))*(-2*x
(1))+2*M*x
(1);
g(2,1)=1+2*M*(-x
(1)^2+x
(2));
end%题8
(1)
functionp=f_pfun(x)
%罚函数
a=-x
(1)^2+x
(2);
b=x
(1);
ifa>=0&b>=0
p=0;
elseifa>=0&b<0
p=b^2;
elseifa<0&b>=0
p=a^2;
else
p=a^2+b^2;
end%题8
(1)
运行结果如下:
x=
1.0e-004*
-0.2058
-0.2058
f=
-2.0576e-005
k=
5
trace=
100.0000-0.0050-0.0017-0.0006-0.0002-0.0001-0.0000
3.0000-0.0050-0.0017-0.0006-0.0002-0.0001-0.0000
从上述运行结果可以得出:
最优解为x=,最小值约为f*=0。
(2)minf(x)=x1^2+x2^2
S.t2x1+x2-2≤0
-x2+1≤0
用matlab求解如下:
主程序与
(1)相同;
注意:
以下三个函数体在运行计算习题四第8
(2)题时使用
%*****************************************************
functiong=f_dfun(x)
globalM;
a=-2*x
(1)-x
(2)+2;
b=x
(2)-1;
ifa>=0&b>=0
g(1,1)=2*x
(1);
g(2,1)=2*x
(2);
elseifa<0&b>=0
g(1,1)=2*x
(1)+2*M*(-2*x
(1)-x
(2)+2)*(-2);
g(2,1)=2*x
(2)+2*M*(-2*x
(1)-x
(2)+2)*(-1);
elseifa>=0&b<0
g(1,1)=2*x
(1);
g(2,1)=2*x
(2)+2*M*(x
(2)-1);
else
g(1,1)=2*x
(1)+2*M*(-2*x
(1)-x
(2)+2)*(-2);
g(2,1)=2*x
(2)+2*M*(-2*x
(1)-x
(2)+2)*(-1)+2*M*(x
(2)-1);
end%题8
(2)
%*****************************************************
functionf=f_fun(x)
globalM;
f=x
(1)^2+x
(2)^2+M*f_pfun(x);%题8
(2)
functionp=f_pfun(x)
a=-2*x
(1)-x
(2)+2;
b=x
(2)-1;
ifa>=0&b>=0
p=0;
elseifa>=0&b<0
p=b^2;
elseifa<0&b>=0
p=a^2;
else
p=a^2+b^2;
end%题8
(2)
运行结果如下:
x=
-0.0000
1.0000
f=
1.0000
k=
6
trace=
100.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000
3.00000.99010.99670.99890.99960.99991.00001.0000
从上述运行结果可以得出:
最优解为x=,最小值约为f*=1。
P232第四章25
用梯度投影法求解下列线性约束优化问题
minf(x)=x1^2+x1*x2+2*x2^2+2*x3^2+2*x2*x3+4*x1+6*x2+12*x3
s.tx1+x2+x3≤6
﹣x1﹣x2+2x3≥2
xi≥0,i=1,2,3
取x1=,
用matlab求解如下:
f=‘x1^2+x1*x2+2*x2^2+2*x3^2+2*x2*x3+4*x1+6*x2+12*x3’;
a=[111;11-2];
b=[6;-2];
l=zeros(3,1);
x0=[113];
[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(f,x0,a,b,[],[],l,[])
运行结果如下:
x=
001
fval=
14
exitflag=
1
output=
iterations:
2
funcCount:
14
stepsize:
1
algorithm:
'medium-scale:
SQP,Quasi-Newton,line-search'
firstorderopt:
0
cgiterations:
[]
message:
[1x144char]
lambda=
lower:
[3x1double]
upper:
[3x1double]
eqlin:
[0x1double]
eqnonlin:
[0x1double]
ineqlin:
[2x1double]
ineqnonlin:
[0x1double]
grad=
4.0000
8.0000
16.0000
hessian=
1.11460.67710.6042
0.67713.36462.4792
0.60422.47923.4583
从上述运行结果可以得出:
最优解为x=,最小值约为f*=14。
P234第四章34
(1)
用乘子法求解
maxf(x)=10x1+4.4x2^2+2x3
s.tx1+4x2+5x3≤32
x1+3x2+2x3≤29
x3^2/2+x2^2≥3
x1≥2,x2≥0,x3≥0
用matlab求解如下:
function[x,minf]=ymh434(l)
formatlong;
symsx1x2x3
f=10*(x1)+4.4*(x2)^2+2*(x3);
h=[-x1-4*x2-5*x3+32,-x1-3*x2-2*x3+29,((x3)^2)/2+(x2)^2-3,x1-2,x2,x3];
x0=[220];
v=[100000];
M=2;
alpha=2;
gama=0.25;
var=[x1x2x3];
eps=1.0e-4;
ifnargin==8
eps=1.0e-4;
end
m1=transpose(x0);
m2=inf;
whilel
FE=0;
u=subs(h,{x1,x2,x3},[m1
(1),m1
(2),m1(3)]);
fori=1:
length(h)
if(v(i)+M*u(i))<=0
FE=FE+((v(i)+M*u(i))^2-(v(i))^2);
else
FE=FE+(v(i))^2;
end
end
SumF=f+(1/(2*M))*FE;
[m2,minf]=minNT(SumF,transpose(m1),var,eps);
Hm2=subs(h,{x1,x2,x3},[m2
(1),m2
(2),m2(3)]);
Hm1=subs(h,{x1,x2,x3},[m1
(1),m1
(2),m1(3)]);
Hx2=Funval(h,var,x2);
Hx1=Funval(h,var,x1);
ifnorm(Hx2)x=x2;
break;
else
ifHx2/Hx1>=gama
M=alpha*M;
x1=x2;
else
v=v-M*transpose(Hx2);
x1=x2;
end
end
end
minf=Funval(f,var,x);
formatshort;
运行结果如下:
x=
19.7502.45
minf=
-202.42
maxf1=
202.42
从上述运行结果可以得出:
最优解为x=,最大值约为f*=202.42,最小值约为f*=-202.42。
P235第四章35
(2)
用序列二次规划法求解
minf(x)=-5x1-5x2-4x3-x1x3-6x4-5x5/(1+x5)-8x6/(1+x6)-10(1-2e^-x7+e^-2x7)
s.tg1(x)=2x4+x5+0.8x6+x7-5=0
g2(x)=x2^2+x3^2+x5^2+x6^2-5=0
g3(x)=x1+x2+x3+x4+x5+x6+x7≤10
g4(x)=x1+x2+x3+x4≤5
g5(x)=x1+x3+x5+x6^2-x7^2-5≤0
xi≥0,i=1,2,…,7
用matlab求解如下:
functionf=objfun35(x)
f=-5*x
(1)-5*x
(2)-4*x(3)-x
(1)*x(3)-6*x(4)-5*x(5)/(1+x(5))-8*x(6)/(1+x(6))-10*(1-2*exp(-x(7))+exp(-2*x(7)));
function[c,ceq]=confun35(x)
ceq=[x
(2)^2+x(3)^2+x(5)^2+x(6)^2-5];
c=[x
(1)+x(3)+x(5)+x(6)^2-x(7)^2-5];
clc
clear
x0=[1,1,1,1,1,1,1];
aeq=[0,0,0,2,1,0.8,1];
beq=5;
a=[1,1,1,1,1,1,1;1,1,1,1,0,0,0];
b=[10,5]';
lb=[0,0,0,0,0,0,0];
ub=[];
[x,fval,exitflag,output,lambda,grad,hessian]=fmincon('objfun35',x0,a,b,aeq,beq,lb,ub,'confun35',options)
[c,ceq]=confun35(x)
运行结果如下:
IterF-countf(x)constraintStep-sizederivativeoptimalityProcedure
08-31.49581Infeasiblestartpoint
117-41.71750.86421-6.31.62
226-43.03090.53581-1.090.924
335-44.51990.83481-1.40.438
444-44.45430.0462110.1260.247
553-44.46360.0043541-0.004560.127
662-44.46970.0053241-0.00070.00717
771-44.46871.631e-00510.0009840.000523
880-44.46876.735e-00913.03e-0067.25e-006
989-44.46874.308e-01311.25e-0094.11e-007
x=3.24180.00001.63420.12400.88961.24022.8702
fval=-44.4687
……
hessian=
0.57380.3772-0.23690.1036-0.1250-0.0757-0.1739
0.37720.68980.47690.0641-0.0816-0.0879-0.0803
-0.23690.47691.45900.06280.06570.1229-0.0055
0.10360.06410.06280.88110.07210.01200.0699
-0.1250-0.08160.06570.07211.7848-0.2129-0.0244
-0.0757-0.08790.12290.0120-0.21291.4472-0.1761
-0.1739-0.0803-0.00550.0699-0.0244-0.17611.0268
c=-5.9342
ceq=4.3077e-013
从上述运行结果可以得出:
最优解为x=,最小值约为f*=-44.4687。