共轭梯度法编程.docx
- 文档编号:9889922
- 上传时间:2023-02-07
- 格式:DOCX
- 页数:20
- 大小:18.45KB
共轭梯度法编程.docx
《共轭梯度法编程.docx》由会员分享,可在线阅读,更多相关《共轭梯度法编程.docx(20页珍藏版)》请在冰豆网上搜索。
共轭梯度法编程
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'*pk a=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。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 共轭 梯度 编程