(完整word版)优化设计Matlab编程作业
(完整word版)优化设计Matlab编程作业
优化设计无约束优化min f(x)= 21x +22x -21x 2x -41x初选x0=[1,1]程序：Step 1: Write an M-file objfun1.m.function f1=objfun1(x)f1=x(1)^2+2*x(2)^2-2*x(1)*x(2)-4*x(1);Step 2: Invoke one of the unconstrained optimization routinesx0=[1,1];>> options = optimset('LargeScale','off');>> [x,fval,exitflag,output] = fminunc(@objfun1,x0,options)运行结果：x =4.0000 2.0000fval =-8.0000exitflag =1output =iterations: 3funcCount: 12stepsize: 1firstorderopt: 2.3842e-007algorithm: 'medium-scale: Quasi-Newton line search'message: [1x85 char]非线性有约束优化1. Min f(x)=321x +22x +21x -32x +5Subject to: 1g (x)=1x +2x +18≤02g (x)=51x -32x -25≤03g (x)=131x -4122x 0≤4g (x)=14≤1x 130≤5g (x)=2≤2x 57≤初选x0=[10,10]Step 1: Write an M-file objfun2.mfunction f2=objfun2(x)f2=3*x(1)^2+x(2)^2+2*x(1)-3*x(2)+5;Step 2: Write an M-file confun1.m for the constraints.function [c,ceq]=confun1(x)% Nonlinear inequality constraintsc=[x(1)+x(2)+18;5*x(1)-3*x(2)-25;13*x(1)-41*x(2)^2;14-x(1);x(1)-130;2-x(2);x(2)-57];% Nonlinear inequality constraintsceq=[];Step 3: Invoke constrained optimization routinex0=[10,10]; % Make a starting guess at the solution>> options = optimset('LargeScale','off');>> [x, fval] = ...fmincon(@objfun2,x0,[],[],[],[],[],[],@confun1,options)运行结果：x =3.6755 -7.0744fval =124.14952. min f （x ）=222154x x +s.t. 0632)(211≤-+=x x x g01)(212≥+=x x x g 初选x0=[1,1]Step 1: Write an M-file objfun3.mfunction f=objfun3(x)f=4*x(1)^2+5*x(2)^2Step 2: Write an M-file confun3.m for the constraints.function [c,ceq]=confun3(x)%Nonlinear inequality constraintsc=[2*x(1)+3*x(2)-6;-x(1)*x(2)-1];% Nonlinear equality constraintsceq[];Step 3: Invoke constrained optimization routinex0=[1,1];% Make a starting guess at the solution>> options = optimset('LargeScale','off');>> [x, fval] = ...fmincon(@objfun,x0,[],[],[],[],[],[],@confun,options)运行结果：Optimization terminated: no feasible solution found. Magnitude of searchdirection less than 2*options.TolX but constraints are not satisfied.x =1 1fval =-13实例：螺栓连接的优化设计图示为一压气机气缸与缸盖连接的示意图。
大连理工优化方法大作业MATLAB编程
function [x,dk,k]=fjqx(x,s) flag=0;a=0;b=0;k=0;d=1;while(flag==0)[p,q]=getpq(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;x=x+d*s;flag=1;endk=k+1;if(p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endend%定义求函数值的函数fun，当输入为x0=（x1，x2）时，输出为f function f=fun(x)f=(x(2)-x(1)^2)^2+(1-x(1))^2;function gf=gfun(x)gf=[-4*x(1)*(x(2)-x(1)^2)+2*(x(1)-1),2*(x(2)-x(1)^2)]; function [p,q]=getpq(x,d,s)p=fun(x)-fun(x+d*s)+0.20*d*gfun(x)*s';q=gfun(x+d*s)*s'-0.60*gfun(x)*s';结果：x=[0,1];s=[-1,1];[x,dk,k]=fjqx(x,s)x =-0.0000 1.0000dk =1.1102e-016k =54function f= fun( X )%所求问题目标函数f=X(1)^2-2*X(1)*X(2)+2*X(2)^2+X(3)^2+ X(4)^2-X(2)*X(3)+2*X(1)+3*X(2)-X(3);endfunction g= gfun( X )%所求问题目标函数梯度g=[2*X(1)-2*X(2)+2,-2*X(1)+4*X(2)-X(3)+3,2*X(3)-X(2)-1,2*X(4)];endfunction [ x,val,k ] = frcg( fun,gfun,x0 )%功能：用FR共轭梯度法求无约束问题最小值%输入：x0是初始点，fun和gfun分别是目标函数和梯度%输出：x、val分别是最优点和最优值，k是迭代次数maxk=5000;%最大迭代次数rho=0.5;sigma=0.4;k=0;eps=10e-6;n=length(x0);while(k<maxk)g=feval(gfun,x0);%计算梯度itern=k-(n+1)*floor(k/(n+1));itern=itern+1;%计算搜索方向if(itern==1)d=-g;elsebeta=(g*g')/(g0*g0');d=-g+beta*d0;gd=g'*d;if(gd>=0.0)d=-g;endendif(norm(g)<eps)break;endm=0;mk=0;while(m<20)if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d) mk=m;break;endm=m+1;endx0=x0+rho^mk*d;val=feval(fun,x0);g0=g;d0=d;k=k+1;endx=x0;val=feval(fun,x0);end结果：>> x0=[0,0,0,0];>> [ x,val,k ] = frcg( 'fun','gfun',x0 )x =-4.0000 -3.0000 -1.0000 0val =-8.0000k =21或者function [x,f,k]=second(x)k=0;dk=dfun(x);g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);while(norm(g1)>=0.02)if(k==3)k=0;g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);else if(k<3)u=((norm(g1))^2)/(norm(g0)^2); s=-g1+u*s;k=k+1;g0=g1;dk=dfun(x);x=x+dk*s;g1=gfun(x);endendf=fun(x);endfunction f=fun(x)f=x(1)^2-2*x(1)*x(2)+2*x(2)^2+x(3)^2+x(4)^2-x(2)*x(3)+2*x(1)+3*x(2)-x(3); function gf=gfun(x)gf=[2*x(1)-2*x(2)+2,-2*x(1)+4*x(2)-x(3)+3,2*x(3)-x(2)-1,2*x(4)];function [p,q]=con(x,d)ss=-gfun(x);p=fun(x)-fun(x+d*ss)+0.2*d*gfun(x)*(ss)';q=gfun(x+d*ss)*(ss)'-0.6*gfun(x)*(ss)';function dk=dfun(x)flag=0;a=0;d=1;while(flag==0)[p,q]=con(x,d);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;flag=1;endif(p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endEnd结果：x=[0,0,0,0];>> [x,f,k]=second(x)x =-4.0147 -3.0132 -1.0090 0 f = -7.9999k = 1function [f,x,k]=third_1(x)k=0;g=gfun(x);while(norm(g)>=0.001)s=-g;dk=dfun(x,s);x=x+dk*s;k=k+1;g=gfun(x);f=fun(x);endfunction f=fun(x)f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)];function [j_1,j_2]=con(x,d,s)j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)'; j_2=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)'; function dk=dfun(x,s)%获取步长flag=0;a=0;d=1;while(flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;flag=1;endif(p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endend结果：x=[0,1];[f,x,k]=third_1(x)f =0.7729x = -0.4196 0.0001k =8（1）程序：function [f,x,k]=third_2(x)k=0;H=inv(ggfun(x));g=gfun(x);while(norm(g)>=0.001)s=(-H*g')';dk=dfun(x,s);x=x+dk*s;k=k+1;g=gfun(x);f=fun(x);endfunction f=fun(x)f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2); function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)]; function ggf=ggfun(x)ggf=[(4*x(1)^2+2)*exp(x(1)^2+x(2)^2),4*x(1)*x(2)*exp(x(1)^2+x(2)^2);4*x(1)*x(2)*exp(x(1)^2+x(2)^2),4+(4*x(2)^2+2)*exp(x(1)^2+x(2)^2)]; function [j_1,j_2]=con(x,d,s)j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)';j_2=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)';function dk=dfun(x,s)% 步长获取flag=0;a=0;d=1;b=10000;while(flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;flag=1;endif(p>=0)&&(q<0)a=d;if 2*d>=(d+b)/2d=(d+b)/2;else d=2*d;endendEnd结果：x=[0,1];[f,x,k]=third_2(x)f =0.7729x = -0.4193 0.0001k =8（2）程序：function [f,x,k]=third_3(x) k=0;X=cell(2);g=cell(2);X{1}=x;H=eye(2);g{1}=gfun(X{1});s=(-H*g{1}')';dk=dfun(X{1},s);X{2}=X{1}+dk*s;g{2}=gfun(X{2});while(norm(g{2})>=0.001)dx=X{2}-X{1};dg=g{2}-g{1};v=dx/(dx*dg')-(H*dg')'/(dg*H*dg'); h1=H*dg'*dg*H/(dg*H*dg');h2=dx'*dx/(dx*dx');h3=dg*H*dg'*v'*v;H=H-h1+h2+h3;k=k+1;X{1}=X{2};g{1}=gfun(X{1});s=(-H*g{1}')';dk=dfun(X{1},s);X{2}=X{1}+dk*s;g{2}=gfun(X{2});norm(g{2});f=fun(x);x=X{2};endfunction f=fun(x)f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)];function ggf=ggfun(x)ggf=[(4*x(1)^2+2)*exp(x(1)^2+x(2)^2),4*x(1)*x(2)*exp(x(1)^2+x(2)^2);4*x(1)*x(2)* exp(x(1)^2+x(2)^2),4+(4*x(2)^2+2)*exp(x(1)^2+x(2)^2);function [p,q]=con(x,d,s)p=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)';q=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)';function dk=dfun(x,s)flag=0;a=0;d=1;b=10000;while(flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0) dk=d;flag=1;endif(p>=0)&&(q<0)a=d;if 2*d>=(d+b)/2d=(d+b)/2;else d=2*d;endendend结果：x=[0,1];[f,x,k]=third_3(x)f =0.7729x = -0.4195 0.0000 k=6function callqpactH=[2 0; 0 2];c=[-2 -5]';Ae=[ ]; be=[ ];Ai=[1 -2; -1 -2; -1 2;1 0;0 1];bi=[-2 -6 -2 0 0]';x0=[0 0]';[x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) function [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) epsilon=1.0e-9; err=1.0e-6;k=0; x=x0; n=length(x); kmax=1.0e3;ne=length(be); ni=length(bi); lamk=zeros(ne+ni,1); index=ones(ni,1);for (i=1:ni)if(Ai(i,:)*x>bi(i)+epsilon), index(i)=0; endendwhile(k<=kmax)Aee=[];if(ne>0), Aee=Ae; endfor(j=1:ni)if(index(j)>0), Aee=[Aee; Ai(j,:)]; end endgk=H*x+c;[m1,n1] = size(Aee);[dk,lamk]=qsubp(H,gk,Aee,zeros(m1,1)); if(norm(dk)<=err)y=0.0;if(length(lamk)>ne)[y,jk]=min(lamk(ne+1:length(lamk))); endif(y>=0)exitflag=0;elseexitflag=1;for(i=1:ni)if(index(i)&(ne+sum(index(1:i)))==jk) index(i)=0; break;endendendk=k+1;elseexitflag=1;alpha=1.0; tm=1.0;for(i=1:ni)if((index(i)==0)&(Ai(i,:)*dk<0)) tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk); if(tm1<tm)tm=tm1; ti=i;endendendalpha=min(alpha,tm);x=x+alpha*dk;if(tm<1), index(ti)=1; end endif(exitflag==0), break; endk=k+1;endoutput.fval=0.5*x'*H*x+c'*x; output.iter=k;function [x,lambda]=qsubp(H,c,Ae,be) ginvH=pinv(H);[m,n]=size(Ae);if(m>0)rb=Ae*ginvH*c + be;lambda=pinv(Ae*ginvH*Ae')*rb;x=ginvH*(Ae'*lambda-c);elsex=-ginvH*c;lambda=0;end结果>>callqpactx =1.40001.7000lambda =0.8000exitflag =output =fval: -6.4500iter: 7function [x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0)%功能: 用乘子法解一般约束问题: min f(x), s.t. h(x)=0, g(x).=0%输入: x0是初始点, fun, dfun分别是目标函数及其梯度；% hf, dhf分别是等式约束（向量）函数及其Jacobi矩阵的转置；% gf, dgf分别是不等式约束（向量）函数及其Jacobi矩阵的转置；%输出: x是近似最优点，mu, lambda分别是相应于等式约束和不等式约束的乘子向量; % output是结构变量, 输出近似极小值f, 迭代次数, 内迭代次数等maxk=500;c=2.0;eta=2.0;theta=0.8;k=0;ink=0;epsilon=0.00001;x=x0;he=feval(hf,x);gi=feval(gf,x);n=length(x);l=length(he);m=length(gi);mu=zeros(l,1);lambda=zeros(m,1);btak=10;btaold=10;while(btak>epsilon&&k<maxk)%调用BFGS算法程序求解无约束子问题[x,ival,ik]=bfgs('mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c);ink=ink+ik;he=feval(hf,x);gi=feval(gf,x);btak=0;for i=1:lbtak=btak+he(i)^2;end%更新乘子向量for i=1:mtemp=min(gi(i),lambda(i)/c);btak=btak+temp^2;endbtak=sqrt(btak);if btak>epsilonif k>=2&&btak>theta*btaoldc=eta*c;endfor i=1:lmu(i)=mu(i)-c*he(i);endfor i=1:mlambda(i)=max(0,lambda(i)-c*gi(i));endk=k+1;btaold=btak;x0=x;endendf=feval(fun,x);output.fval=f;output.iter=k;%增广拉格朗日函数function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c) f=feval(fun,x);he=feval(hf,x);gi=feval(gf,x);l=length(he);m=length(gi);psi=f;s1=0;for i=1:lpsi=psi-he(i)*mu(i);s1=s1+he(i)^2;endpsi=psi+0.5*c*s1;s2=0;for i=1:ms3=max(0,lambda(i)-c*gi(i));s2=s2+s3^2-lambda(i)^2;endpsi=psi+s2/(2*c);%不等式约束函数文件g1.mfunction gi=g1(x)gi=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;%目标函数的梯度文件df1.mfunction g=df1(x)g=[4, -2*x(2)]';%等式约束（向量）函数的Jacobi矩阵（转置）文件dh1.m function dhe=dh1(x)dhe=[-2*x(1), -2*x(2)]'%不等式约束（向量）函数的Jacobi矩阵（转置）文件dg1.m function dgi=dg1(x)dgi=[10-2*x(1), 10-2*x(2)]';function [x,val,k]=bfgs(fun,gfun,x0,varargin)maxk=500;rho=0.55;sigma=0.4;epsilon=0.00001;k=0;n=length(x0);Bk=eye(n);while(k<maxk)gk=feval(gfun,x0,varargin{:});if(norm(gk)<epsilon)break;enddk=-Bk\gk;m=0;mk=0;while(m<20)newf=feval(fun,x0+rho^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<oldf+sigma*rho^m*gk'*dk)mk=m;break;endm=m+1;endx=x0+rho^mk*dk;sk=x-x0;yk=feval(gfun,x,varargin{:})-gk;if(yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;x0=x;endval=feval(fun,x0,varargin{:});结果x=[2 2]';[x,mu,lambda,output]=multphr('fun','hf','gf1','df','dh','dg',x0) x =1.00134.8987mu =0.7701lambda =0.9434output =fval: -31.9923iter: 4f=[3,1,1];A=[2,1,1;1,-1,-1];b=[2;-1];lb=[0,0,0];x=linprog(f,A,b,zeros(3),[0,0,0]',lb)结果：Optimization terminated.x =0.00000.50000.5000。
(完整word版)MATLAB程序设计教程(第二版)课后答案
MATLAB第二版课后答案unit3—8 unit3实验指导1、 n=input（'请输入一个三位数：'）;a=fix（n/100)；b=fix（（n-a*100)/10）；c=n-a＊100-b＊10；d=c＊100+b*10+a2（1）n=input('请输入成绩’）；switch ncase num2cell（90：100)p='A'；case num2cell（80:89)p='B'；case num2cell(70：79）p=’C';case num2cell（60:69)p='D'；otherwisep='E';endprice=p(2）n=input(’请输入成绩');if n〉=90&n〈=100p='A’；elseif n>=80&n<=89p='B';elseif n〉=70&n<=79p=’C’;elseif n>=60＆n<=69p='D'；elsep='E'；endprice=p(3）tryn;catchprice='erroe'end3n=［1,5,56，4,3,476，45，6，3，76，45,6,4,3，6，4,23,76,908，6]; a=n(1);b=n(1)；for m=2:20if n（m）>aa=n（m）；elseif n（m)<bb=n（m）;endendmax=amin=b法2n=[1,5，56,4，3，476,45，6,3,76,45，6，4，3,6，4，23，76,908,6］;min=min(n)max=max(n）4b=［—3.0:0.1:3.0];for n=1：61a=b(n）；y(n)=（exp（0.3＊a)-exp(—0。
(完整word版)matlab经典习题及解答
第1章 MATLAB概论1.1与其他计算机语言相比较，MATLAB语言突出的特点是什么？MATLAB具有功能强大、使用方便、输入简捷、库函数丰富、开放性强等特点.1.2 MATLAB系统由那些部分组成?MATLAB系统主要由开发环境、MATLAB数学函数库、MATLAB语言、图形功能和应用程序接口五个部分组成。
1.4 MATLAB操作桌面有几个窗口?如何使某个窗口脱离桌面成为独立窗口?又如何将脱离出去的窗口重新放置到桌面上?在MATLAB操作桌面上有五个窗口,在每个窗口的右上角有两个小按钮，一个是关闭窗口的Close按钮，一个是可以使窗口成为独立窗口的Undock按钮，点击Undock按钮就可以使该窗口脱离桌面成为独立窗口,在独立窗口的view菜单中选择Dock ……菜单项就可以将独立的窗口重新防止的桌面上。
1.5 如何启动M文件编辑/调试器？在操作桌面上选择“建立新文件”或“打开文件"操作时，M文件编辑/调试器将被启动.在命令窗口中键入edit命令时也可以启动M文件编辑/调试器.1.6 存储在工作空间中的数组能编辑吗？如何操作？存储在工作空间的数组可以通过数组编辑器进行编辑：在工作空间浏览器中双击要编辑的数组名打开数组编辑器，再选中要修改的数据单元，输入修改内容即可。
1。
7 命令历史窗口除了可以观察前面键入的命令外,还有什么用途?命令历史窗口除了用于查询以前键入的命令外，还可以直接执行命令历史窗口中选定的内容、将选定的内容拷贝到剪贴板中、将选定内容直接拷贝到M文件中。
1。
8 如何设置当前目录和搜索路径,在当前目录上的文件和在搜索路径上的文件有什么区别？当前目录可以在当前目录浏览器窗口左上方的输入栏中设置，搜索路径可以通过选择操作桌面的file菜单中的Set Path 菜单项来完成。
在没有特别说明的情况下，只有当前目录和搜索路径上的函数和文件能够被MATLAB运行和调用，如果在当前目录上有与搜索路径上相同文件名的文件时则优先执行当前目录上的文件，如果没有特别说明，数据文件将存储在当前目录上。
Matlab优化设计作业答案
长江大学机械工程学院机械优化设计大作业班级2012年5月31-361、⎩⎨⎧=+-=++⋅+-++=22422min 321321321232221x x x x x x t s x x x x x x f解： function f=fun1(x)f=x(1)^2+2*x(2)^2+x(3)^2-2*x(1)*x(2)+x(3)clearclcx0=[1;1;1];A=[];b=[];Aeq=[1 1 1;2 -1 1];beq=[4;2];lb=[];ub=[];[x,fval,exitflag]=fmincon('fun1',x0,A,b,Aeq,beq,lb,ub)x =1.90911.95450.1364fval =3.9773exitflag =12、221212min ()2130f x x x x s t x =+-+⋅-≤解： function f=fun2(x)f=x(1)^2+x(2)^2-2*x(1)+1clearclcx0=[1;1];A=[0 -1];b=-3;Aeq=[];beq=[];lb=[];ub=[];[x,fval,exitflag]=fmincon('fun2',x0,A,b,Aeq,beq,lb,ub)x =13fval =9exitflag =13、 2212min ()(2)f X x x =-+112222312..()0()0()10s t g X x g X x g X x x =-≤=-≤=-+-≤解：function f=fun3(x)f=(x(1)-2)^2+x(2)^2function [g,ceq]=mycon3(x)g=[-x(1);-x(2);-x(1)^2+x(2)^2-1];ceq=[];clearclcx0=[0;0];A=[];b=[];Aeq=[];beq=[];lb=[];ub=[];[x,fval]=fmincon('fun3',x0,A,b,Aeq,beq,lb,ub,'mycon3')x =2.0000fval =2.2204e-0164、⎩⎨⎧≤--≤+--⋅++++=01005.1)12424(min 21212122122211x x x x x x t s x x x x x e f x解： function f=fun4(x)f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);function [g,ceq]=mycon4(x)g(1)=1.5+x(1)*x(2)-x(1)-x(2);g(2)=-x(1)*x(2)-10;ceq=[];clearclcx0=[0;0];A=[];b=[];Aeq=[];beq=[];lb=[];ub=[];[x,fval]=fmincon('fun4',x0,A,b,Aeq,beq,lb,ub,'mycon4')x =1.1825-1.7398fval =3.06085 喜糖问题：需要购买甲乙两种喜糖，喜糖甲10元/斤，喜糖乙20元/斤。
(完整word版)Matlab求解优化问题
预备知识：M 文件简介在MATLAB 中，用户可以利用Edtior （编辑器）建立M 文件，然后在命令窗口中的“>>”提示符下键入M 文件的主文件名，回车执行.MATLAB 的M 文件有两类：命令文件和函数文件。
将原本要在MATLAB 环境下直接输入的语句，放在一个以 .m 为后缀的文件中，这一文件就称为命令文件；函数文件由五部分组成：函数定义行、H1行、函数帮助文本、函数体、注释，MATLAB 的内部函数都是由函数文件定义的。
1.11 优化(最值、数学规划)在数学上，优化问题包括最值问题和数学规划问题等，后者又包括线性规划、整数规划（含0-1规划）、二次规划等.在MATLAB 中，求解最值问题的命令主要有：fminbnd (f,x1,x2) 求一元函数f 在区间[x1,x2]上的最小值点[x,fval]=fminbnd(f,x1,x2) 求一元函数f 在区间[x1,x2]上的最小值点和最小值 fminsearch (’f’,x0) 求多元函数f 在点x0附近的最小值点[x,fval]=fminsearch(’f’,x0) 求多元函数f 在点x0附近的最小值点和最小值例1.11.1 求函数23)(2++=x x x f 在区间]5,5[-上的最小值点和最小值. >> [x,fval]=fminbnd('x^2+3*x+2',-5,5) x =-1.5000 fval =-0.2500例1.11.2 求函数21212122),(x x x x x x f ++=在点)1,1(附近的最小值点和最小值. >> [x,fval]= fminsearch('x(1)*x(2)+2/x(1)+2/x(2)',[1 1]) x =1.2599 1.2599 fval =4.7622在MATLAB 中，求解数学规划问题的命令主要有：（1）线性规划⎪⎪⎩⎪⎪⎨⎧≤≤=⋅≤=ub x lb beq x Aeq bAx t s x c z T ..min命令：[x,fval]=linprog (c,A,b,Aeq,beq,lb,ub)在上述命令中，当某些参数空缺时，可用[]代替或省略，下同。
(完整word版)matlab(7,4)汉明码和(7,4)循环码的编程设计
二、创新实验设计创新实验一：（7,4）汉明码的编码与译码实现1、实验目的实现（7,4）汉明码的编码与译码，通过这次实验不但加深了对汉明码编码和译码原理了解，而且对线性分组码有所了解。
2、实验原理线性分组码的构造方法比较简单、理论较为成熟，应用比较广泛。
汉明码是一种能够纠正一个错码的效率比较高的线性分组码，下面以（7,4）码为例就汉明码的编码与译码分别进行介绍：（1）编码原理一般来说，若汉明码长为n ，信息位数为k ，则监督位数r=n-k 。
若希望用r 个监督位构造出r 个监督关系式来指示一位错码的n 种可能位置，则要求21r n -≥或211rk r -≥++ (1)设汉明码（n,k ）中k=4，为了纠正一位错码，由式（1）可知，要求监督位数r ≥3。
若取r=3,则n=k+r=7。
这样就构成了（7,4）码。
用6543210a a a a a a a 来表示这7个码元，用123s s s 的值表示3个监督关系式中的校正子，则123s s s 的值与错误码元位置的对应关系可以规定如表1所列。
表2.1 校正子和错码位置的关系则由表1可得监督关系式：16542s a a a a =⊕⊕⊕()226531s a a a a =⊕⊕⊕()3 36430s a a a a =⊕⊕⊕()4 在发送端编码时，信息位6543a a a a 的值决定于输入信号，因此它们是随机的。
监督位2a 、1a 、a 应根据信息位的取值按监督关系来确定，为使所编的码中无错码，则123,,S S S 等于0，即65426531643000(5)0a a a a a a a a a a a a ⊕⊕⊕=⎧⎪⊕⊕⊕=⎨⎪⊕⊕⊕=⎩方程组（5）可等效成如下矩阵形式6543210111010001101010010110010a a a a a a a ⎡⎤⎢⎥⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎢⎥⎣⎦（6）式(6)可简化为0T T HA =，H 为监督矩阵，则由式(6)可得到监督矩阵11101001101010=[P I ] (7)1011001r H ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦因为生成矩阵'=[I Q]=[I ]k k G P ，所以由(7)得生成矩阵G 如下：[]k 10001110100110[']00101010001011k G I Q I P ⎡⎤⎢⎥⎢⎥===⎢⎥⎢⎥⎣⎦然后利用信息位和生成矩阵G 相乘产生整个码组，即有[][]65432106543=(8)A a a a a a a a a a a a G=其中A 为整个码组矩阵，6543a a a a 是信息位。
最优化方法matlab作业
实用最优化方法——matlab编程作业初值为[-1;1]其中g0、g1分别为不同x值下得导数，f0、f1为函数值MATLAB程序：x0=[-1;1];s0=[1;1];c1=0.1;c2=0.5;a=0;b=inf;d=1;n=0;x1=x0+d*s0;g0=[-400*(x0(2)-x0(1)^2)*x0(1)-2*(1-x0(1));200*(x0(2)-x0(1) ^2)];g1=[-400*(x1(2)-x1(1)^2)*x1(1)-2*(1-x1(1));200*(x1(2)-x1(1) ^2)];f1=100*(x1(2)-x1(1)^2)^2+(1-x1(1))^2;f0=100*(x0(2)-x0(1)^2)^2+(1-x0(1))^2;while((f0-f1<-c1*d*g0'*s0)||(g1'*s0<c2*g0'*s0))if ((f0-f1)<(-c1*d*g0'*s0))b=d;d=(d+a)/2;x1=x0+d*s0;g0=[-400*(x0(2)-x0(1)^2)*x0(1)-2*(1-x0(1));200*(x0 (2)-x0(1)^2)];g1=[-400*(x1(2)-x1(1)^2)*x1(1)-2*(1-x1(1));200*(x1 (2)-x1(1)^2)];f1=100*(x1(2)-x1(1)^2)^2+(1-x1(1))^2;f0=100*(x0(2)-x0(1)^2)^2+(1-x0(1))^2;elseif (((g1')*s0)<(c2*(g0')*s0))a=d;if(2*d<=(d+b)/2)d=2*d;elsed=(d+b)/2;endx1=x0+d*s0;g0=[-400*(x0(2)-x0(1)^2)*x0(1)-2*(1-x0(1));200*(x0(2) -x0(1)^2)];g1=[-400*(x1(2)-x1(1)^2)*x1(1)-2*(1-x1(1));200*(x1(2 )-x1(1)^2)];f1=100*(x1(2)-x1(1)^2)^2+(1-x1(1))^2;f0=100*(x0(2)-x0(1)^2)^2+(1-x0(1))^2;endx1df1=100*(x1(2)-x1(1)^2)^2+(1-x1(1))^2 计算结果：最优点：x1 = -0.99611.0039步长： d = 0.0039最优解：f1 = 3.9981目标函数：function f = fun2( x )f=x(1)^2-2*x(1)*x(2)+2*x(2)^2+x(3)^2-x(2)*x(3)+2*x(1)+3*x(2 )-x(3);end目标函数梯度：function g = gfun2( x )g=[2 -2 0;-2 4 -1;0 -1 2]*x+[2;3;-1];end源代码：（以初值为为（0；0；0））x0=[0;0;0]; %初始值eps=1.0e-5; %精度g0=gfun2(x0);s0=-g0;n=0;syms d1;while norm(g0)>epsif n<3g=gfun2(x0+d1*s0);d= double(solve(s0'*g));x1=x0+d*s0;g1=gfun2(x1);if norm(g1)<epsn=n+1;x0=x1;breakelses0=-g1+(norm(g1)^2/norm(g0)^2)*s0;x0=x1;g0=g1;endelseif n==3x0=x1;g0=gfun2(x0);s0=-g0;n=0;endn=n+1;x0nfun2(x0)计算结果：最优点：x0 = -4 -3 -1 迭代次数： n = 3 最优值： ans = -8（1）最速下降法：目标函数：function f= fun3_1(x )f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);end目标函数的梯度：function g= gfun3_1(x)g=[1+2*x(1)*exp(x(1)^2+x(2)^2);4*x(2)+2*x(2)*exp(x(1)^2+x(2 )^2)];end源代码（初值为（1；0））：x0=[1;0];%初始值eps=1.0e-5;%精度n=0;g0=gfun3(x0);syms d1;while norm(g0)>=epss0=-g0;g=gfun3(x0+d1*s0);d= double(solve(s0'*g));x1=x0+d*s0;g1=gfun3(x1);if( norm(g1)<eps)n=n+1;x0=x1;break;elsex0=x1;g0=gfun3(x0);endn=n+1;endf0=fun3(x0)x0n计算结果：最优值： f0 = 0.7729 最优点：x0 = -0.4194 0迭代次数：n = 1（2）牛顿法目标函数：function f= fun3_1(x )f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);end目标函数梯度：function g= gfun3_1(x)g=[1+2*x(1)*exp(x(1)^2+x(2)^2);4*x(2)+2*x(2)*exp(x(1)^2+x(2 )^2)];end目标函数的Hesse阵：function g2 = g2fun3(x)g2=[2*exp(x(1)^2+x(2)^2)+4*x(1)^2*exp(x(1)^2+x(2)^2),4*x(1) *x(2)*exp(x(1)^2+x(2)^2)4*x(1)*x(2)*exp(x(1)^2+x(2)^2),4+2*exp(x(1)^2+x(2)^2)+4*x(2 )^2*exp(x(1)^2+x(2)^2)];end源代码（初值为（1；0））：x0=[1;0];%初始值eps=1.0e-5;%精度n=0;g0=gfun3(x0);g20=g2fun3(x0);while norm(g0)>=epsd=-g20\g0;x1=x0+d;g1=gfun3(x1);if( norm(g1)<eps)n=n+1;x0=x1;break;elsex0=x1;g0=gfun3(x0);endn=n+1;endf0=fun3(x0)x0n计算结果：最优值：f0 = 0.7729最优点：x0 = -0.4194迭代次数：n = 63（3）利用BGFS法：目标函数：function f= fun3_1(x )f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);end目标函数梯度function g= gfun3_1(x)g=[1+2*x(1)*exp(x(1)^2+x(2)^2);4*x(2)+2*x(2)*exp(x(1)^2+x(2 )^2)];end源代码（初值为（1；0））：x0=[1;0];%初始值eps=1.0e-5;%精度n=0;g0=gfun3(x0);syms d1;h0=eye(2);while norm(g0)>=epss0=-h0*g0;g=gfun3(x0+d1*s0);d= double(solve(s0'*g));x1=x0+d*s0;g1=gfun3(x1);if( norm(g1)<eps)n=n+1;x0=x1;break;elseh0=h0-(h0*(g1-g0)*(g1-g0)'*h0)/((g1-g0)'*h0*(g1-g0))+...((x1-x0)*(x1-x0)')/((x1-x0)'*(g1-g0))+...((g1-g0)'*h0*(g1-g0))*((x1-x0)*(x1-x0)');x0=x1;g0=gfun3(x0);endn=n+1;endf0=fun3(x0)x0n计算结果：最优值：f0 = 0.7729最优点：x0 = -0.4194迭代次数：n = 1题四：求解子程序：function [x,lambda]=qsubp(H,c,Ae,be) ginvH=pinv(H);[m,n]=size(Ae);if(m>0)rb=Ae*ginvH*c + be;lambda=pinv(Ae*ginvH*Ae')*rb;x=ginvH*(Ae'*lambda-c);elsex=-ginvH*c;lambda=0;endend源代码：H=[2 0;0 2];%目标函数的hesse阵c=[-2,-4];Ae=[0 0];be=[0;0];Ai=[1/2 0;-1 3];bi=[1;2];x0=[-1;0];%初始值内部点eps=1.0e-9; %µ±ax-b=epsʱµ±×öax-b=0´¦Àíerr=1.0e-6;k=0;x0(1)=x0(1)+x0(2);x=x0;n=length(x);max=1.0e3;ne=length(be);ni=length(bi);lamk=zeros(ne+ni,1);index=ones(ni,1);for i=(1:ni)if(bi(i)-Ai(i,:)*x>eps)index(i)=0;。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
优化设计
无约束优化
min f(x)= 21x +22x -21x 2x -41x
初选x0=[1,1]
程序：
Step 1: Write an M-file objfun1.m.
function f1=objfun1(x)
f1=x(1)^2+2*x(2)^2-2*x(1)*x(2)-4*x(1);
Step 2: Invoke one of the unconstrained optimization routines
x0=[1,1];
>> options = optimset('LargeScale','off');
>> [x,fval,exitflag,output] = fminunc(@objfun1,x0,options)
运行结果：
x =
4.0000 2.0000
fval =
-8.0000
exitflag =
1
output =
iterations: 3
funcCount: 12
stepsize: 1
firstorderopt: 2.3842e-007
algorithm: 'medium-scale: Quasi-Newton line search'
message: [1x85 char]
非线性有约束优化
1. Min f(x)=321x +2
2x +21x -32x +5
Subject to: 1g (x)=1x +2x +18≤0
2g (x)=51x -32x -25≤0
3g (x)=131x -412
2x 0≤
4g (x)=14≤1x 130≤
5g (x)=2≤2x 57≤
初选x0=[10,10]
Step 1: Write an M-file objfun2.m
function f2=objfun2(x)
f2=3*x(1)^2+x(2)^2+2*x(1)-3*x(2)+5;
Step 2: Write an M-file confun1.m for the constraints.
function [c,ceq]=confun1(x)
% Nonlinear inequality constraints
c=[x(1)+x(2)+18;
5*x(1)-3*x(2)-25;
13*x(1)-41*x(2)^2;
14-x(1);
x(1)-130;
2-x(2);
x(2)-57];
% Nonlinear inequality constraints
ceq=[];
Step 3: Invoke constrained optimization routine
x0=[10,10]; % Make a starting guess at the solution
>> options = optimset('LargeScale','off');
>> [x, fval] = ...
fmincon(@objfun2,x0,[],[],[],[],[],[],@confun1,options)
运行结果：
x =
3.6755 -7.0744
fval =
124.1495
2. min f （x ）=2
22154x x +
s.t. 0632)(211≤-+=x x x g
01)(212≥+=x x x g 初选x0=[1,1]
Step 1: Write an M-file objfun3.m
function f=objfun3(x)
f=4*x(1)^2+5*x(2)^2
Step 2: Write an M-file confun3.m for the constraints.
function [c,ceq]=confun3(x)
%Nonlinear inequality constraints
c=[2*x(1)+3*x(2)-6;
-x(1)*x(2)-1];
% Nonlinear equality constraints
ceq[];
Step 3: Invoke constrained optimization routine
x0=[1,1];% Make a starting guess at the solution
>> options = optimset('LargeScale','off');
>> [x, fval] = ...
fmincon(@objfun,x0,[],[],[],[],[],[],@confun,options)
运行结果：
Optimization terminated: no feasible solution found. Magnitude of search
direction less than 2*options.TolX but constraints are not satisfied.
x =
1 1
fval =
-13
实例：螺栓连接的优化设计
图示为一压气机气缸与缸盖连接的示意图。
已知D1=400mm,D2=240mm ，缸内工作压力p=8.5Mpa ，螺栓材料为45Cr ，抗拉强度Mpa 1000b =σ，屈服强度Mpa 320s =σ，拉压疲劳极限Mpa 3301-=σ，许用疲劳安全系数[7.1]S a =，取残余预紧力F F 6.1'
'=，采用铜皮石棉密封垫片，螺栓相对刚度8.0=c K 。
从安全、可靠、经济的角度来选择螺栓的个数n 和螺栓的直径d 。
解：
1. 目标函数
取螺栓组连接经济成本n C 最小为目标。
当螺栓的长度、材料和加工条件一定时，螺栓的总成本与n ，d 值成正比，故本问题优化设计的目标函数为 min 21)(x x nd C x f n ===
由此可见，设计变量为螺栓个数n 和直径d 为
n x [= 1[]x d T = T x ]2
2．约束条件
（1）强度约束条件：螺栓在脉动载荷下工作，因此螺栓组连接须满足疲劳强度条件
()()()[]a
a a S K K S ≥++-+=-min min 122σσψσψσσσσσ 其中a σ为应力幅值;min σ为最小应力；σK 为疲劳极限综合影响系数，取σK =4.4；σψ为应力折算系数，取σψ=0.23. 气缸最大载荷p D P 224
1π=
螺栓最大工作载荷n P F F +=,,0，210max 4
1d F πσ+ 螺栓最小工作载荷n p F F K F F /8.000,
-=-=σ, 21,
min 41D F πσ= 螺栓应力幅值2min
max σσσ-=a
对于普通螺纹，小径d d 85.01=，于是疲劳强度约束条件为
()[]01≤-=a a S S X g
（2）密封约束条件：考虑密封安全，螺栓间距应小于8d ，故密封约束条件为
()08400821
1
2≤-=-=x x d n D X g ππ （3）扳手工作空间约束条件：考虑扳手工作空间，螺距间距应大于2d ，故扳手工作空间约束条件为()0400221
213≤-=-
=x x n D d X g ππ （4）非约束条件 (),011≤-=x X g ()025≤-=x X g
螺栓连接的优化数学模型
综上所述，本问题的数学模型可表达如下
设计变量：[]T
x x x 21,= 目标函数：()21m in x x x f =
约束条件：()5,4,3,2,1,0..=≤u x g t s u
M 文件
function f=stud_obj(x)
f=x(1)*x(2);
global p Ksigam psai sigam_1
p=8.5*le6;
Ksigam=4.4;psai=0.23;sigam_1=330;D2=240;
p=1/4*pi*D2^2*le-6*p;
function [c,ceq]=stud_conl(x)
global p Ksigam psai sigam_1
p,x
F0=(1.6*p+p)/x(1);
d1=0.85*x(2);
A=1/4*pi*d1^2;
F1=F0-0.8*p/x(1);
sigma_max=F0/A;
sigma_min=F1/A;
sigma_1=(sigma_max-sigma_min)/2;
Sa=(2*sigma_1+(Ksigma-psai)*sigma_min)/((Ksigma+psai)*(2*sigma_a+sigm a_min));
c=[1.7-Sa;
400*pi/x(1)-8*x(2);
-400*pi/x(1)+2*x(2);
16-x(1);
-x(2)];
ceq=[];
x0=[7,20];
[x,feval]=fmincon(@stud_obj,x0,[],[],[],[],[],[],@stud_conl)
运行结果：
x =
16.00 28.84
fval =
461.39
根据实际问题的意义取整、标准化：n=16，d=30，经验证n和d的取值满足约束条件。