代數(shù)方程 最優(yōu)化方程

代數(shù)方程最優(yōu)化方程第一頁，共六十一頁，2022年，8月28日7.1代數(shù)方程的求解

7.1.1代數(shù)方程的圖解法一元方程的圖解法例：>>ezplot('exp(-3*t)…*sin(4*t+2)+4*exp…(-0.5*t)*cos(2*t)-…0.5',[05])>>holdon,>>line([0,5],[0,0])%同時(shí)繪制橫軸第二頁，共六十一頁，2022年，8月28日驗(yàn)證：>>symst;t=3.52028;vpa(exp(-3*t)*sin(4*t+2)+…4*exp(-0.5*t)*cos(2*t)-0.5)ans=第三頁，共六十一頁，2022年，8月28日二元方程的圖解法例：>>ezplot('x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y)')%第一個(gè)方程曲線>>holdon%保護(hù)當(dāng)前坐標(biāo)系>>ezplot(‘x^2*…cos(x+y^2)+…y^2*exp(x+y)')第四頁，共六十一頁，2022年，8月28日方程的圖解法僅適用于一元、二元方程的求根問題。第五頁，共六十一頁，2022年，8月28日7.1.2多項(xiàng)式型方程的準(zhǔn)解析解法例：>>ezplot('x^2+y^2-1');holdon%繪制第一方程的曲線>>ezplot(‘0.75*x^3-y+0.9’)%繪制第二方程為關(guān)于x的6次多項(xiàng)式方程應(yīng)有6對(duì)根。圖解法只能顯示求解方程的實(shí)根。第六頁，共六十一頁，2022年，8月28日一般多項(xiàng)式方程的根可為實(shí)數(shù)，也可為復(fù)數(shù)?？捎肕ATLAB符號(hào)工具箱中的solve()函數(shù)。最簡調(diào)用格式：

S=solve(eqn1,eqn2,…,eqnn)

（返回一個(gè)結(jié)構(gòu)題型變量S，如S.x表示方程的根。）直接得出根：（變量返回到MATLAB工作空間）

[x,…]=solve(eqn1,eqn2,…,eqnn)同上，并指定變量

[x,…]=solve(eqn1,eqn2,…,eqnn，’x,…’)第七頁，共六十一頁，2022年，8月28日例：>>symsxy;>>[x,y]=solve('x^2+y^2-1=0','75*x^3/100-y+9/10=0')x=y=第八頁，共六十一頁，2022年，8月28日驗(yàn)證>>[eval('x.^2+y.^2-1')eval('75*x.^3/100-y+9/10')]ans=[0,0][0,0][0,0][-.1e-31,0][.5e-30+.1e-30*i,0][.5e-30-.1e-30*i,0]

由于方程階次可能太高，不存在解析解。然而，可利用MATLAB的符號(hào)工具箱得出原始問題的高精度數(shù)值解，故稱之為準(zhǔn)解析解。第九頁，共六十一頁，2022年，8月28日例：>>[x,y,z]=solve('x+3*y^3+2*z^2=1/2','x^2+3*y+z^3=2','x^3+2*z+2*y^2=2/4');>>size(x)ans=271>>x(22),y(22),z(22)ans=ans=.37264668450644375527750811296627e-1ans=第十頁，共六十一頁，2022年，8月28日驗(yàn)證：>>err=[x+3*y.^3+2*z.^2-1/2,x.^2+3*y+z.^3-2,x.^3+2*z+2*y.^2-2/4];>>norm(double(eval(err)))ans=1.4998e-027多項(xiàng)式乘積形式也可，如把第三個(gè)方程替換一下。>>[x,y,z]=solve('x+3*y^3+2*z^2=1/2','x^2+3*y+z^3=2','x^3+2*z*y^2=2/4');>>err=[x+3*y.^3+2*z.^2-1/2,x.^2+3*y+z.^3-2,x.^3+2*z.*y.^2-2/4];>>norm(double(eval(err)))%將解代入求誤差ans=5.4882e-028第十一頁，共六十一頁，2022年，8月28日例：求解（含變量倒數(shù)）>>symsxy;>>[x,y]=solve('x^2/2+x+3/2+2/y+5/(2*y^2)+3/x^3=0',...'y/2+3/(2*x)+1/x^4+5*y^4','x,y');>>size(x)ans=261>>err=[x.^2/2+x+3/2+2./y+5./(2*y.^2)+3./x.^3,y/2+3./(2*x)+1./x.^4+5*y.^4];％驗(yàn)證>>norm(double(eval(err)))ans=8.9625e-030第十二頁，共六十一頁，2022年，8月28日例：求解（帶參數(shù)方程）>>symsabxy;>>[x,y]=solve('x^2+a*x^2+6*b+3*y^2=0','y=a+(x+3)','x,y')x=[1/2/(4+a)*(-6*a-18+2*(-21*a^2-45*a-27-24*b-6*a*b-3*a^3)^(1/2))][1/2/(4+a)*(-6*a-18-2*(-21*a^2-45*a-27-24*b-6*a*b-3*a^3)^(1/2))]y=[a+1/2/(4+a)*(-6*a-18+2*(-21*a^2-45*a-27-24*b-6*a*b-3*a^3)^(1/2))+3][a+1/2/(4+a)*(-6*a-18-2*(-21*a^2-45*a-27-24*b-6*a*b-3*a^3)^(1/2))+3]第十三頁，共六十一頁，2022年，8月28日7.1.3一般非線性方程數(shù)值解格式：最簡求解語句

x=fsolve(Fun,x0)一般求解語句

[x,f,flag,out]=fsolve(Fun,x0,opt,p1,p2,…)若返回的flag大于0，則表示求解成功，否則求解出現(xiàn)問題，opt求解控制參數(shù)，結(jié)構(gòu)體數(shù)據(jù)。獲得默認(rèn)的常用變量opt=optimset;用這兩種方法修改參數(shù)（解誤差控制量）opt.Tolx=1e-10;或set(opt.‘Tolx’,1e-10)第十四頁，共六十一頁，2022年，8月28日例：自編函數(shù)：functionq=my2deq(p)q=[p(1)*p(1)+p(2)*p(2)-1;0.75*p(1)^3-p(2)+0.9];>>OPT=optimset;OPT.LargeScale='off';>>[x,Y,c,d]=fsolve('my2deq',[1;2],OPT)Optimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.x=0.35700.9341Y=1.0e-009*0.12150.0964第十五頁，共六十一頁，2022年，8月28日c=1d=iterations:7funcCount:21algorithm:'trust-regiondogleg'firstorderopt:1.3061e-010％解回代的精度調(diào)用inline()函數(shù)：>>f=inline('[p(1)*p(1)+p(2)*p(2)-1;0.75*p(1)^3-p(2)+0.9]','p');>>[x,Y]=fsolve(f,[1;2],OPT);%結(jié)果和上述完全一致，從略。Optimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.若改變初始值x0=[-1,0]T第十六頁，共六十一頁，2022年，8月28日>>[x,Y,c,d]=fsolve(f,[-1,0]',OPT);x,Y,kk=d.funcCountOptimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.x=-0.98170.1904Y=1.0e-010*0.5619-0.4500kk=15初值改變有可能得出另外一組解。故初值的選擇對(duì)解的影響很大，在某些初值下甚至無法搜索到方程的解。第十七頁，共六十一頁，2022年，8月28日例：用solve()函數(shù)求近似解析解>>symst;solve(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5)ans=％不允許手工選擇初值，只能獲得這樣的一個(gè)解?？上扔糜脠D解法選取初值，再調(diào)用fsolve()函數(shù)數(shù)值計(jì)算第十八頁，共六十一頁，2022年，8月28日>>formatlong>>y=inline('exp(-3*t).*sin(4*t+2)+4*exp(-0.5*t).*cos(2*t)-0.5','t');>>ff=optimset;ff.Display='iter';[t,f]=fsolve(y,3.5203,ff)NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius121.8634e-0095.16e-0051243.67694e-0193.61071e-0057.25e-0101Optimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.t=3.52026389294877f=-6.063776702980306e-010第十九頁，共六十一頁，2022年，8月28日重新設(shè)置相關(guān)的控制變量，提高精度。>>ff=optimset;ff.TolX=1e-16;ff.TolFun=1e-30;>>ff.Display='iter';[t,f]=fsolve(y,3.5203,ff)NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius121.8634e-0095.16e-0051243.67694e-0193.61071e-0057.25e-01013605.07218e-01001Optimizationterminatedsuccessfully:First-orderoptimalityislessthanoptions.TolFun.t=3.52026389244155f=0第二十頁，共六十一頁，2022年，8月28日7.2無約束最優(yōu)化問題求解

7.2.1解析解法和圖解法第二十一頁，共六十一頁，2022年，8月28日例：>>symst;y=exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5;>>y1=diff(y,t);%求取一階導(dǎo)函數(shù)>>ezplot(y1,[0,4])%繪制出選定區(qū)間內(nèi)一階導(dǎo)函數(shù)曲線第二十二頁，共六十一頁，2022年，8月28日>>t0=solve(y1)%求出一階導(dǎo)數(shù)等于零的點(diǎn)t0=>>y2=diff(y1);b=subs(y2,t,t0)%并驗(yàn)證二階導(dǎo)數(shù)為正b=>>t=0:0.4:4;y=exp(-3*t).*sin(4*t+2)+4*exp(-0.5*t).*cos(2*t)-0.5;>>plot(t,y)第二十三頁，共六十一頁，2022年，8月28日7.2.2基于MATLAB的數(shù)值解法第二十四頁，共六十一頁，2022年，8月28日例：>>f=inline('(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2))','x');>>x0=[0;0];ff=optimset;ff.Display='iter';>>x=fminsearch(f,x0,ff)IterationFunc-countminf(x)Procedure13-0.000499937initial24-0.000499937reflect

72137-0.641424contractoutsideOptimizationterminatedsuccessfully:x=0.6111-0.3056第二十五頁，共六十一頁，2022年，8月28日

>>x=fminunc(f,[0;.0],ff)DirectionalIterationFunc-countf(x)Step-sizederivative12-2e-0080.001-429-0.5846690.3043530.343316-0.6413970.9503220.00191422-0.6414240.984028-1.45e-008x=0.6110-0.3055比較可知fminunc()函數(shù)效率高于fminsearch()函數(shù)，故若安裝了最優(yōu)化工具箱則應(yīng)調(diào)用fminunc()函數(shù)。第二十六頁，共六十一頁，2022年，8月28日7.2.3全局最優(yōu)解與局部最優(yōu)解例：>>f=inline('exp(-2*t).*cos(10*t)+exp(-3*(t+2)).*sin(2*t)','t');%目標(biāo)函數(shù)>>t0=1;[t1,f1]=fminsearch(f,t0);[t1f1]ans=0.9228-0.1547>>t0=0.1;[t2,f2]=fminsearch(f,t0);[t2f2]ans=0.2945-0.5436第二十七頁，共六十一頁，2022年，8月28日>>symst;y=exp(-2*t)*cos(10*t)+exp(-3*(t+2))*sin(2*t);>>ezplot(y,[0,2.5]);set(gca,‘Ylim’,[-0.6,1])％在t[0,2.5]內(nèi)的曲線>>ezplot(y,[-0.5,2.5]);set(gca,‘Ylim’,[-2,1.2])％在[-0.5,2.5]曲線>>t0=-0.2;[t3,f3]=fminsearch(f,t0);[t3f3]ans=-0.3340-1.9163第二十八頁，共六十一頁，2022年，8月28日7.2.4利用梯度求解最優(yōu)化問題例：>>[x,y]=meshgrid…(0.5:0.01:1.5);…z=100*(y.^2-x).^2…+(1-x).^2;>>contour3(x,y,z,100),set(gca,'zlim',[0,310])%測試算法的函數(shù)第二十九頁，共六十一頁，2022年，8月28日>>f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2','x');>>ff=optimset;ff.TolX=1e-10;ff.TolFun=1e-20;ff.Display='iter';>>x=fminunc(f,[0;0],ff)Warning:Gradientmustbeprovidedfortrust-regionmethod;usingline-searchmethodinstead.DirectionalIterationFunc-countf(x)Step-sizederivative1210.5-4442023.01749e-0123.40397e-009-1.77e-013x=1.000001736959721.00000347608069第三十頁，共六十一頁，2022年，8月28日％求梯度向量(比較引入梯度的作用，但梯度的計(jì)算也費(fèi)時(shí)間）>>symsx1x2;f=100*(x2-x1^2)^2+(1-x1)^2;>>J=jacobian(f,[x1,x2])J=[-400*(x2-x1^2)*x1-2+2*x1,200*x2-200*x1^2]function[y,Gy]=c6fun3(x)％目標(biāo)函數(shù)編寫：y=100*(x(2)-x(1)^2)^2+(1-x(1))^2;%目標(biāo)函數(shù)Gy=[-400*(x(2)-x(1)^2)*x(1)-2+2*x(1);200*x(2)-200*x(1)^2];%梯度>>ff.GradObj='on';x=fminunc('c6fun3',[0;0],ff)NormofFirst-orderIterationf(x)stepoptimalityCG-iterations196.38977e-0292.12977e-0121.6e-0141x=0.999999999999990.99999999999998第三十一頁，共六十一頁，2022年，8月28日7.3有約束最優(yōu)化問題的計(jì)算機(jī)求解

6.3.1約束條件與可行解區(qū)域有約束最優(yōu)化問題的一般描述：對(duì)于一般的一元問題和二元問題，可用圖解法直接得出問題的最優(yōu)解。第三十二頁，共六十一頁，2022年，8月28日例：用圖解方法求解：>>[x1,x2]=meshgrid(-3:.1:3);%生成網(wǎng)格型矩陣>>z=-x1.^2-x2;%計(jì)算出矩陣上各點(diǎn)的高度>>i=find(x1.^2+x2.^2>9);z(i)=NaN;%找出x1^2+x2^2>9的坐標(biāo)，并置函數(shù)值為NaN>>i=find(x1+x2>1);z(i)=NaN;%找出x1+x2>1的坐標(biāo),置為NaN>>surf(x1,x2,z);shadinginterp;>>max(z(:)),view(0,90)ans=3第三十三頁，共六十一頁，2022年，8月28日7.3.2線性規(guī)劃問題的計(jì)算機(jī)求解第三十四頁，共六十一頁，2022年，8月28日例：求解>>f=-[21431]';A=[02142;345-1-1];>>B=[54;62];Ae=[];Be=[];xm=[0,0,3.32,0.678,2.57];>>ff=optimset;ff.LargeScale='off';%不使用大規(guī)模問題求解>>ff.TolX=1e-15;ff.TolFun=1e-20;ff.TolCon=1e-20;ff.Display='iter';>>[x,f_opt,key,c]=linprog(f,A,B,Ae,Be,xm,[],[],ff)第三十五頁，共六十一頁，2022年，8月28日Optimizationterminatedsuccessfully.x=19.78500.00003.320011.38502.5700f_opt=-89.5750key=1％求解成功c=iterations:5algorithm:'medium-scale:activeset'cgiterations:[]第三十六頁，共六十一頁，2022年，8月28日例：求解>>f=[-3/4,150,-1/50,6]';Aeq=[];Beq=[];>>A=[1/4,-60,-1/50,9;1/2,-90,-1/50,3];B=[0;0];>>xm=[-5;-5;-5;-5];xM=[Inf;Inf;1;Inf];ff=optimset;>>ff.TolX=1e-15;ff.TolFun=1e-20;TolCon=1e-20;ff.Display='iter';>>[x,f_opt,key,c]=linprog(f,A,B,Aeq,Beq,xm,xM,[0;0;0;0],ff)第三十七頁，共六十一頁，2022年，8月28日Residuals:PrimalDualUpperDualityTotalInfeasInfeasBoundsGapRelA*x-bA'*y+z-w-f{x}+s-ubx'*z+s'*wError-------------------------------------------------------------Iter0:9.39e+0031.43e+0021.94e+0026.03e+0042.77e+001Iter1:6.38e-0121.21e+0010.00e+0003.50e+0031.78e+000

Iter10:0.00e+0006.15e-0260.00e+0001.70e-0244.10e-028Optimizationterminatedsuccessfully.x=-5.0000-0.19471.0000-5.0000f_opt=-55.4700key=1c=iterations:10cgiterations:0algorithm:'lipsol'第三十八頁，共六十一頁，2022年，8月28日7.3.3二次型規(guī)劃的求解第三十九頁，共六十一頁，2022年，8月28日例：求解>>f=[-2,-4,-6,-8];H=diag([2,2,2,2]);>>OPT=optimset;OPT.LargeScale='off';%不使用大規(guī)模問題算法>>A=[1,1,1,1;3,3,2,1];B=[5;10];Aeq=[];Beq=[];LB=zeros(4,1);>>[x,f_opt]=quadprog(H,f,A,B,Aeq,Beq,LB,[],[],OPT)Optimizationterminatedsuccessfully.x=0.00000.66671.66672.6667f_opt=-23.6667％（解的目標(biāo)值應(yīng)為30＋（-23.6667）＝6.3333)第四十頁，共六十一頁，2022年，8月28日6.3.4一般非線性規(guī)劃問題的求解第四十一頁，共六十一頁，2022年，8月28日例：求解目標(biāo)函數(shù)：functiony=opt_fun1(x)y=1000-x(1)*x(1)-2*x(2)*x(2)-x(3)*x(3)-x(1)*x(2)-x(1)*x(3);非線性約束函數(shù)：function[c,ceq]=opt_con1(x)ceq=[x(1)*x(1)+x(2)*x(2)+x(3)*x(3)-25;8*x(1)+14*x(2)+7*x(3)-56];c=[];第四十二頁，共六十一頁，2022年，8月28日>>ff=optimset;ff.LargeScale='off';ff.Display='iter';>>ff.TolFun=1e-30;ff.TolX=1e-15;ff.TolCon=1e-20;>>x0=[1;1;1];xm=[0;0;0];xM=[];A=[];B=[];Aeq=[];Beq=[];>>[x,f_opt,c,d]=fmincon('opt_fun1',x0,A,B,Aeq,Beq,xm,xM,'opt_con1',ff)；>>x,f_opt,kk=d.funcCountx=3.51210.21703.5522f_opt=961.7152kk=%目標(biāo)函數(shù)調(diào)用的次數(shù)108第四十三頁，共六十一頁，2022年，8月28日簡化非線約束函數(shù)function[c,ceq]=opt_con2(x)ceq=x(1)*x(1)+x(2)*x(2)+x(3)*x(3)-25;c=[];求解：>>x0=[1;1;1];Aeq=[8,14,7];Beq=56;>>[x,f_opt,c,d]=fmincon('opt_fun1',x0,A,B,Aeq,Beq,xm,xM,'opt_con2',ff);maxDirectionalFirst-orderIterF-countf(x)constraintStep-sizederivativeoptimalityProcedure111955.33622.90.25-29518.3infeasible21116961.71501-6.3e-0156.97e-005HessianmodifiedtwiceOptimizationterminatedsuccessfully:Searchdirectionlessthan2*options.TolXandmaximumconstraintviolationislessthanoptions.TolConActiveConstraints:12>>x,f_opt,kk=d.funcCount%從略（計(jì)算結(jié)果同上一樣）第四十四頁，共六十一頁，2022年，8月28日例：利用梯度求解梯度函數(shù)：>>symsx1x2x3;f=1000-x1*x1-2*x2*x2-x3*x3-x1*x2-x1*x3;>>J=jacobian(f,[x1,x2,x3])J=[-2*x1-x2-x3,-4*x2-x1,-2*x3-x1]改寫目標(biāo)函數(shù)：function[y,Gy]=opt_fun2(x)y=1000-x(1)*x(1)-2*x(2)*x(2)-x(3)*x(3)-x(1)*x(2)-x(1)*x(3);Gy=-[2*x(1)+x(2)+x(3);4*x(2)+x(1);2*x(3)+x(1)];第四十五頁，共六十一頁，2022年，8月28日>>x0=[1;1;1];xm=[0;0;0];xM=[];A=[];B=[];Aeq=[];Beq=[];>>ff=optimset;ff.GradObj=‘on’;％若知道梯度函數(shù)ff.Display='iter';ff.LargeScale='off';>>ff.TolFun=1e-30;ff.TolX=1e-15;ff.TolCon=1e-20;>>[x,f_opt,c,d]=fmincon('opt_fun2',x0,A,B,Aeq,Beq,xm,xM,'opt_con1',ff);>>x,f_opt,kk=d.funcCountx=3.51210.21703.5522f_opt=961.7152kk=95第四十六頁，共六十一頁，2022年，8月28日7.4整數(shù)規(guī)劃問題的計(jì)算機(jī)求解7.4.1整數(shù)線性規(guī)劃問題的求解A、B線性等式和不等式約束，約束式子由ctype變量控制，intlist為整數(shù)約束標(biāo)示，how＝0表示結(jié)果最優(yōu)，2為無可行解，其余失敗。第四十七頁，共六十一頁，2022年，8月28日例：>>f=-[21431]';A=[02142;345-1-1];intlist=ones(5,1);>>B=[54;62];ctype=[-1;-1];xm=[0,0,3.32,0.678,2.57];xM=inf*ones(5,1);>>[res,b]=ipslv_mex(f,A,B,intlist,xM,xm,ctype)%因?yàn)榉祷氐腷=0，表示優(yōu)化成功res=1904105b=0第四十八頁，共六十一頁，2022年，8月28日>>[x1,x2,x3,x4,x5]=ndgrid(1:20,0:20,4:20,1:20,3:20);>>i=find((2*x2+x3+4*x4+2*x5<=54)&(3*x1+4*x2+5*x3-x4-x5<=62));>>f=-2*x1(i)-x2(i)-4*x3(i)-3*x4(i)-x5(i);[fmin,ii]=sort(f);>>index=i(ii(1));x=[x1(index),x2(index),x3(index),x4(index),x5(index)]x=1904105％當(dāng)把20換為30，一般計(jì)算機(jī)配置下實(shí)現(xiàn)不了，故求解整數(shù)規(guī)劃時(shí)不適合采用窮舉算法。第四十九頁，共六十一頁，2022年，8月28日次最優(yōu)解>>fmin(1:10)'ans=-89-88-88-88-88-88-88-88-87-87>>in=i(ii(1:4));x=[x1(in),x2(in),x3(in),x4(in),x5(in),fmin(1:4)]x=1904105-891804113-881705104-881506104-88>>intlist=[1,0,0,1,1];％混合整數(shù)規(guī)劃>>[res,b]=ipslv_mex(f,A,B,intlist,xM,xm,ctype)res=19.000003.800011.00003.0000b=0第五十頁，共六十一頁，2022年，8月28日7.4.2一般非線性整數(shù)規(guī)劃問題與求解第五十一頁，共六十一頁，2022年，8月28日err字符串為空，則返回最優(yōu)解。該函數(shù)尚有不完全之處，解往往不是精確整數(shù)，可用下面語句將其化成整數(shù)。第五十二頁，共六十一頁，2022年，8月28日例：functionf=c6miopt(x)f=-[21431]*x;>>A=[02142;345-1-1];intlist=ones(5,1);Aeq=[];Beq=[];>>B=[54;62];ctype=[-1;-1];>>xm=[0,0,4,1,3]';xM=20000*ones(5,1);x0=xm;>>[errmsg,f,X]=bnb20('c6miopt',x0,intlist,xm,xM,A,B,Aeq,Beq);X=X'X=19.000004.000010.00005.0000第五十三頁，共六十一頁，2022年，8月28日>>intlist=[1,0,0,1,1]';>>xm=[0,