最優(yōu)化方法實(shí)驗(yàn)報(bào)告

最優(yōu)化方法實(shí)驗(yàn)報(bào)告optimizationmethodExperimentReport學(xué)生所在學(xué)院：理學(xué)院學(xué)生所在班級：信息1學(xué)生姓名：教務(wù)處20014年5月最優(yōu)化方法實(shí)驗(yàn)報(bào)告書說明:1.下面程序在MATLABR2012a中均能正常運(yùn)行。程序之間有關(guān)聯(lián)。實(shí)驗(yàn)一熟悉MATLAB根本功能〔２學(xué)時〕實(shí)驗(yàn)的目的和要求：在本次實(shí)驗(yàn)中，通過親臨使用MATLAB，對該軟件做一全面了解并掌握重點(diǎn)內(nèi)容。實(shí)驗(yàn)內(nèi)容：全面了解MATLAB系統(tǒng)實(shí)驗(yàn)常用工具的具體操作和功能學(xué)習(xí)建議：本次實(shí)驗(yàn)在全面了解軟件系統(tǒng)根底之上，學(xué)習(xí)和熟悉一些MATLAB的根底用途，重點(diǎn)掌握優(yōu)化工具箱函數(shù)選用的內(nèi)容。重點(diǎn)和難點(diǎn)：優(yōu)化工具箱函數(shù)選用。利用Matlab的優(yōu)化工具箱，可以求解線性規(guī)劃、非線性規(guī)劃和多目標(biāo)規(guī)劃問題。具體而言，包括線性、非線性最小化，最大最小化，二次規(guī)劃，半無限問題，線性、非線性方程〔組〕的求解，線性、非線性的最小二乘問題。另外，該工具箱還提供了線性、非線性最小化，方程求解，曲線擬合，二次規(guī)劃等問題中大型課題的求解方法，為優(yōu)化方法在工程中的實(shí)際應(yīng)用提供了更方便快捷的途徑。優(yōu)化工具箱中的函數(shù)優(yōu)化工具箱中的函數(shù)包括下面幾類：1．最小化函數(shù)表1最小化函數(shù)表函

數(shù)描

述fgoalattain多目標(biāo)到達(dá)問題fminbnd有邊界的標(biāo)量非線性最小化fmincon有約束的非線性最小化fminimax最大最小化fminsearch,fminunc無約束非線性最小化fseminf半無限問題linprog線性課題quadprog二次課題2．方程求解函數(shù)表2方程求解函數(shù)表函

數(shù)描

述/線性方程求解fsolve非線性方程求解fzero標(biāo)量非線性方程求解3．最小二乘〔曲線擬合〕函數(shù)表3最小二乘函數(shù)表函

數(shù)描

述/線性最小二乘lsqlin有約束線性最小二乘lsqcurvefit非線性曲線擬合lsqnonlin非線性最小二乘lsqnonneg非負(fù)線性最小二乘4．實(shí)用函數(shù)表4實(shí)用函數(shù)表函

數(shù)描

述optimset設(shè)置參數(shù)optimget

5．大型方法的演示函數(shù)表5大型方法的演示函數(shù)表函

數(shù)描

述circustent馬戲團(tuán)帳篷問題—二次課題molecule用無約束非線性最小化進(jìn)行分子組成求解optdeblur用有邊界線性最小二乘法進(jìn)行圖形處理6．中型方法的演示函數(shù)表6中型方法的演示函數(shù)表函

數(shù)描

述bandemo香蕉函數(shù)的最小化dfildemo過濾器設(shè)計(jì)的有限精度goaldemo目標(biāo)到達(dá)舉例optdemo演示過程菜單tutdemo教程演示下面以我們最常用的線性規(guī)劃模型求解函數(shù)linprog作為典型對優(yōu)化工具箱進(jìn)行簡單的介紹。linprog函數(shù)功能：求解線性規(guī)劃問題。在命令窗口，鍵入doclinprog,得到下列圖〔該圖為幫助窗口〕數(shù)學(xué)模型：其中f,x,b,beq,lb和ub為向量，A和Aeq為矩陣。語法：x=linprog(f,A,b,Aeq,beq)x=linprog(f,A,b,Aeq,beq,lb,ub)x=linprog(f,A,b,Aeq,beq,lb,ub,x0)x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options)[x,fval]=linprog(...)[x,fval,exitflag]=linprog(...)[x,fval,exitflag,output]=linprog(...)[x,fval,exitflag,output,lambda]=linprog(...)描述：x=linprog(f,A,b)求解問題minf'*x，約束條件為A*x<=b。x=linprog(f,A,b,Aeq,beq)求解上面的問題，但增加等式約束，即Aeq*x=beq。假設(shè)沒有不等式存在，那么令A(yù)=[]、b=[]。x=linprog(f,A,b,Aeq,beq,lb,ub)定義設(shè)計(jì)變量x的下界lb和上界ub，使得x始終在該范圍內(nèi)。假設(shè)沒有等式約束，令A(yù)eq=[]、beq=[]。x=linprog(f,A,b,Aeq,beq,lb,ub,x0)設(shè)置初值為x0。該選項(xiàng)只適用于中型問題，缺省時大型算法將忽略初值。x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options)用options指定的優(yōu)化參數(shù)進(jìn)行最小化。[x,fval]=linprog(...)返回解x處的目標(biāo)函數(shù)值fval。[x,lambda,exitflag]=linprog(...)返回exitflag值，描述函數(shù)計(jì)算的退出條件。[x,lambda,exitflag,output]=linprog(...)返回包含優(yōu)化信息的輸出變量output。[x,fval,exitflag,output,lambda]=linprog(...)將解x處的拉格朗日乘子返回到lambda參數(shù)中。變量：lambda參數(shù)lambda參數(shù)是解x處的拉格朗日乘子。它有以下一些屬性：

lambda.lower–lambda的下界。

lambda.upper–lambda的上界。

lambda.ineqlin–lambda的線性不等式。

lambda.eqlin–lambda的線性等式。其它參數(shù)意義同前。算法：大型優(yōu)化算法

大型優(yōu)化算法采用的是LIPSOL法，該法在進(jìn)行迭代計(jì)算之前首先要進(jìn)行一系列的預(yù)處理。中型優(yōu)化算法

linprog函數(shù)使用的是投影法，就象quadprog函數(shù)的算法一樣。linprog函數(shù)使用的是一種活動集方法，是線性規(guī)劃中單純形法的變種。它通過求解另一個線性規(guī)劃問題來找到初始可行解。診斷：大型優(yōu)化問題

算法的第一步涉及到一些約束條件的預(yù)處理問題。有些問題可能導(dǎo)致linprog函數(shù)退出，并顯示不可行的信息。在本例中，exitflag參數(shù)將被設(shè)為負(fù)值以表示優(yōu)化失敗。假設(shè)Aeq參數(shù)中某行的所有元素都為零，但Beq參數(shù)中對應(yīng)的元素不為零，那么顯示以下退出信息：Exitingduetoinfeasibility:anallzerorowintheconstraintmatrixdoesnothaveazeroincorrespondingrighthandsizeentry.假設(shè)x的某一個元素沒在界內(nèi)，那么給出以下退出信息：

Exitingduetoinfeasibility:objectivef'*xisunboundedbelow.假設(shè)Aeq參數(shù)的某一行中只有一個非零值，那么x中的相關(guān)值稱為奇異變量。這里，x中該成分的值可以用Aeq和Beq算得。假設(shè)算得的值與另一個約束條件相矛盾，那么給出以下退出信息：Exitingduetoinfeasibility:Singletonvariablesinequalityconstraintsarenotfeasible.假設(shè)奇異變量可以求解但其解超出上界或下界，那么給出以下退出信息：Exitingduetoinfeasibility:singletonvariablesintheequalityconstraintsarenotwithinbounds.應(yīng)用實(shí)例這是matlab幫助窗口里給出的一個例子：Findxthatminimizesf(x)=–5x1–4x2–6x3,subjecttox1–x2+x3≤203x1+2x2+4x3≤423x1+2x2≤300≤x1,0≤x2,0≤x3.First,enterthecoefficientsf=[-5;-4;-6];A=[1-11324320];b=[20;42;30];lb=zeros(3,1);Next,callalinearprogrammingroutine.[x,fval,exitflag,output,lambda]=linprog(f,A,b,[],[],lb);ExaminethesolutionandLagrangemultipliers:x,lambda.ineqlin,lambda.lowerx=0.000015.00003.0000ans=0.00001.50000.5000ans=1.00000.00000.0000下面在來用linprog解我們最優(yōu)化考試的題：Minf(x)=-3x1+x2+x3;S.T.x1–2*x2+x3<=11-4*x1+x2+2*x3–x4=3-2*x1+x3=1X1,x2,x3,x4>=0;在matlabcommandwindow中鍵入以下指令：f=[-3;1;1];>>A=[1-21;4-1-2];>>b=[11;-3];>>Aeq=[-201];>>beq=1;>>lb=zeros(3,1);>>[x,fval,exitflag]=linprog(f,A,b,Aeq,beq,lb)Optimizationterminated.x=4.00001.00009.0000fval=-2.0000exitflag=1實(shí)驗(yàn)二一維搜索方法的MATLAB實(shí)現(xiàn)〔２學(xué)時〕實(shí)驗(yàn)的目的和要求：通過本次實(shí)驗(yàn)應(yīng)使學(xué)生掌握如何使用MATLAB軟件進(jìn)行一維搜索，并學(xué)會對具體問題進(jìn)行分析。實(shí)驗(yàn)內(nèi)容：0.618法的MATLAB實(shí)現(xiàn)Fibonacci法的MATLAB實(shí)現(xiàn)學(xué)習(xí)建議：本次實(shí)驗(yàn)是學(xué)生初次使用MATLAB進(jìn)行優(yōu)化問題的實(shí)驗(yàn)，本次實(shí)驗(yàn)就是要通過對一些具體問題的分析學(xué)會軟件的操作并加深對理論知識的理解。重點(diǎn)和難點(diǎn)：具體問題的步長因子確實(shí)定，理解、掌握精度與效率的關(guān)系。實(shí)驗(yàn)內(nèi)容：0.618法和Fibonacci法都是分割方法，其根本思想是通過取試探點(diǎn)和進(jìn)行函數(shù)值的比擬，使包含極小點(diǎn)的搜索區(qū)間不斷縮短，當(dāng)區(qū)間長度縮短到一定程度時，區(qū)間上每個點(diǎn)的函數(shù)值均接近極小值，從而各點(diǎn)可以看作為極小點(diǎn)的近似。這類方法僅需計(jì)算函數(shù)值，不涉及導(dǎo)數(shù)，又稱直接法。他們用途很廣，尤其適用于非光滑及導(dǎo)數(shù)表達(dá)式復(fù)雜或?qū)懖怀龅那樾?。注意，這些方法要求所考慮區(qū)間上的目標(biāo)函數(shù)是單峰函數(shù)，如果這個條件不滿足，我們可以把所考慮的區(qū)間分成假設(shè)干個小區(qū)間，在每個區(qū)間上的函數(shù)式單峰的。這樣，我們在每個小區(qū)間上求極小點(diǎn)，然后選取其中的最小點(diǎn)。一0.618法1.0.618法方法原理:0.618法的根本思想是通過取試探點(diǎn)使包含極小點(diǎn)的區(qū)間(不確定區(qū)間)不斷縮短,當(dāng)區(qū)間長度小到一定程度時,區(qū)間上各點(diǎn)的函數(shù)值均接近極小值,因此任意一點(diǎn)都可作為極小點(diǎn)的近似.0.618法計(jì)算試探點(diǎn)的公式:2.0.618法的算法步驟:=1\*GB3①置初始區(qū)間及精度要求,計(jì)算試探點(diǎn)和,計(jì)算函數(shù)值和.計(jì)算公式是令.=2\*GB3②假設(shè),那么停止計(jì)算.否那么,當(dāng)時,轉(zhuǎn)步驟=3\*GB3③;當(dāng)時,轉(zhuǎn)步驟=4\*GB3④.=3\*GB3③置,,,,計(jì)算函數(shù)值,轉(zhuǎn)步驟=5\*GB3⑤.=4\*GB3④置,,,,計(jì)算函數(shù)值,轉(zhuǎn)步驟=5\*GB3⑤.=5\*GB3⑤置,返回步驟=2\*GB3②.MATLAB實(shí)現(xiàn):3.代碼及數(shù)值算例:程序源代碼:function[X,FMIN,K]=find0618(f,a1,b1,e)%[X,FMIN,K]=find0618(f,a1,b1,e)0.618法一維搜索%f目標(biāo)函數(shù)%a1,b1初始區(qū)間%e精度要求%X極小點(diǎn)%FMIN極小值%K迭代次數(shù)%2014張超a=a1;b=b1;k=1;r=a+0.328*(b-a);u=a+0.618*(b-a);while1iff(r)>f(u)if(b-r)<=eu;break;elsea=r;b=b;r=u;u=a+0.618*(b-a);endelseif(u-a)<=er;break;elsea=a;b=u;u=r;r=a+0.382*(b-a);endk=k+1;endX=(r+u)/2;FMIN=double(f(X));K=k;end數(shù)值算例:Minf(x)=2*x*x–x–1;初始區(qū)間,精度e<=0.16.鍵入命令并輸出結(jié)果：symsxf(x)=2*x^2-x-1;a1=-1;b1=1;e=0.16;[X,FMIN,K]=find0618(f,a1,b1,e)X=0.2258FMIN=-1.1238K=6二Fibonacci法1.Fibonacci法根本原理和步驟思想：搜索區(qū)間長度縮短率采用Fibonacci數(shù)123581321345589……MATLAB實(shí)現(xiàn):2.代碼及數(shù)值算例:程序源代碼:function[X,Fmin,K]=fibonacci(f,a0,b0,e)%fibonacci()Fibonacci法求極小值%X極值點(diǎn)%Fmin極小值%K需要用到第K個Fibonacci數(shù)%a0,b0初始搜索區(qū)間%e精度%張超編寫于2014/04/01a=a0;b=b0;F=[11];i=1;whileF(i)<=(b-a)/eF(i+2)=F(i)+F(i+1);i=i+1;endm=i;r=a+F(m-2)/F(m)*(b-a);u=a+F(m-1)/F(m)*(b-a);fork=1:m-3iff(r)<f(u)a=a;b=u;u=r;r=a+F(m-k-2)/F(m-k)*(b-a);elsea=r;b=b;r=u;u=a+F(m-k-1)/F(m-k)*(b-a);endendX=(r+u)/2;Fmin=double(f(X));K=m;end數(shù)值算例:Minf(x)=x*x–x+2;初始區(qū)間,精度e<=0.08.容易驗(yàn)證，在此區(qū)間上的函數(shù)為嚴(yán)格凸函數(shù)。為了進(jìn)行比擬我們給出其精確解：t*=0.5,f(t*)=1.75。鍵入命令并輸出結(jié)果：symsxf(x)=x^2-x+2;a1=-1;b1=3;e=0.08;[X,FMIN,K]=fibonacci(f,a1,b1,e)X=0.5273FMIN=1.7507K=10實(shí)驗(yàn)三無約束最優(yōu)化方法的MATLAB實(shí)現(xiàn)〔２學(xué)時〕實(shí)驗(yàn)的目的和要求：通過本次實(shí)驗(yàn)使學(xué)生進(jìn)一步熟悉掌握使用MATLAB軟件，并能利用該軟件進(jìn)行無約束最優(yōu)化方法的計(jì)算。實(shí)驗(yàn)內(nèi)容：１、最速下降法的MATLAB實(shí)現(xiàn)２、牛頓法的MATLAB實(shí)現(xiàn)３、共軛梯度法的MATLAB實(shí)現(xiàn)學(xué)習(xí)建議：本次實(shí)驗(yàn)就是要通過對一些具體問題的分析進(jìn)一步熟悉軟件的操作并加深對理論知識的理解。重點(diǎn)和難點(diǎn)：通過同一個具體問題用不同的方法解決的比擬，加深理解恰中選用優(yōu)化問題解決方法的重要性。一最速下降法1.最速下降法根本原理和步驟思想：尋求最速下降方向即負(fù)梯度方向MATLAB實(shí)現(xiàn):2.代碼及數(shù)值算例:程序源代碼:function[X,FMIN,K]=zuisuxiajiang(f,x,x0,e)%[X,FMIN,N]=zuisuxiajiang()法求解無約束問題%X極小點(diǎn)%FMIN極小值%K迭代次數(shù)%f問題函數(shù)%x變量%x0初始點(diǎn)%e終止誤差%張超編寫于2014/04/15count=0;td=jacobian(f,x)';whilenorm(subs(td,x,x0))>eP=-subs(td,x,x0);symsry=x0+r*P;ft(r)=subs(f,x,y);[r0]=fibonacci(ft,0,100,0.01);x0=x0+r0*P;count=count+1;endX=x0;FMIN=subs(f,x,x0);K=count;end二牛頓法1.牛頓法根本原理和步驟思想：在第k次迭代的迭代點(diǎn)x(k)鄰域內(nèi)，用一個二次函數(shù)〔如二階泰勒多項(xiàng)式〕去近似代替原目標(biāo)函數(shù)f(x)，然后求出該二次函數(shù)的極小點(diǎn)作為對原目標(biāo)函數(shù)求優(yōu)的下一個迭代點(diǎn)，依次類推，通過屢次重復(fù)迭代，使迭代點(diǎn)逐步逼近原目標(biāo)函數(shù)的極小點(diǎn)。設(shè)f(x)二次連續(xù)可微，在點(diǎn)x(k〕處的Hesse矩陣正定。MATLAB實(shí)現(xiàn):2.代碼及數(shù)值算例:程序源代碼:function[X,FMIN,K]=ysNewton(f,x,x0,e)%[X,FMIN,N]=ysNewton()原始牛頓法求解無約束問題%X極小點(diǎn)%FMIN極小值%K迭代次數(shù)%f問題函數(shù)%x變量%x0初始點(diǎn)%e終止誤差%張超編寫于2014/04/15count=0;td=jacobian(f,x)';H=jacobian(td',x);whilenorm(subs(td,x,x0))>eP=-subs(H,x,x0)^(-1)*subs(td,x,x0);x0=x0+P;count=count+1;endX=x0;FMIN=subs(f,x,x0);K=count;end牛頓法對于二次正定函數(shù)只需做一次迭代就得到最優(yōu)解。特別在極小點(diǎn)附近，收斂性很好速度也很快。但牛頓法也有缺點(diǎn)，它要求初始點(diǎn)離最優(yōu)解不遠(yuǎn)，假設(shè)初始點(diǎn)選的離最優(yōu)解太遠(yuǎn)時，牛頓法并不能保證其收斂，甚至也不是下降方向。為了克服牛頓法的缺點(diǎn)，人們保存了從牛頓法中選取牛頓方向作為搜索方向，摒棄其步長恒取1的作法，而用一維搜索確定最優(yōu)步長來構(gòu)造算法。程序源代碼:function[X,FMIN,K]=xzNewton(f,x,x0,e)%[X,FMIN,N]=xzNewton()帶步長牛頓法求解無約束問題%X極小點(diǎn)%FMIN極小值%K迭代次數(shù)%f問題函數(shù)%x變量%x0初始點(diǎn)%e終止誤差%張超編寫于2014/04/15count=0;td=jacobian(f,x)';H=jacobian(td',x);whilenorm(subs(td,x,x0))>eP=-subs(H,x,x0)^(-1)*subs(td,x,x0);symsry=x0+r*P;ft(r)=subs(f,x,y);[r0]=fibonacci(ft,0,100,0.01);x0=x0+r0*P;count=count+1;endX=x0;FMIN=subs(f,x,x0);K=count;end三共軛梯度法1.共軛梯度法根本原理和步驟思想：將共軛性和最速下降方向相結(jié)合，利用迭代點(diǎn)處的梯度方向構(gòu)造一組共軛方向，并沿此方向進(jìn)行搜索，求出函數(shù)的極小點(diǎn)。MATLAB實(shí)現(xiàn):2.代碼及數(shù)值算例:程序源代碼:function[X,FMIN,K]=gongetidu(f,x,x0,e)%[X,FMIN,N]=gongetidu()共軛梯度法求解無約束問題%X極小點(diǎn)%FMIN極小值%K迭代次數(shù)%f問題函數(shù)%x變量%x0初始點(diǎn)%e終止誤差%張超編寫于2014/04/15count=1;td=jacobian(f,x)';H=jacobian(td',x);ifnorm(subs(td,x,x0))>eP=-subs(td,x,x0);r0=-subs(td,x,x0)'*P/(P'*H*P);x0=x0+r0*P;elsex0;endwhilenorm(double(subs(td,x,x0)))>eb0=subs(td,x,x0)'*subs(td,x,x0)/(P'*P);P=-subs(td,x,x0)+b0*P;r0=-subs(td,x,x0)'*P/(P'*H*P);x0=x0+r0*P;count=count+1;endX=x0;FMIN=subs(f,x,x0);K=count;end四一個算例分別用上述三中方法計(jì)算下題,并比擬各算法.Minf(x)=(x1-2)^2+(x1–2*x2)^2初始點(diǎn)x0=(0,3)T允許誤差e=0.1鍵入命令并輸出結(jié)果：symsx1x2>>f=(x1-2)^2+(x1-2*x2)^2;>>x=[x1;x2];>>x0=[0;3];>>e=0.1;[X,FMIN,N]=zuisuxiajiang(f,x,x0,e)X=1.97630.9818FMIN=7.2076e-04N=10>>[X,FMIN,N]=ysNewton(f,x,x0,e)X=21FMIN=0N=1[X,FMIN,N]=gongetidu(f,x,x0,e)X=21FMIN=0N=2由上述結(jié)果我們發(fā)現(xiàn)：對于二次正定函數(shù)newton法只需一次迭代就得到正確結(jié)果，共軛梯度法只需進(jìn)行兩次〔因?yàn)槟繕?biāo)函數(shù)是二元函數(shù)〕迭代就得出正確結(jié)果。但最速下降法卻迭代了10次，雖然一維搜索存在誤差，但實(shí)際上最速下降法也需迭代屢次。實(shí)驗(yàn)四約束最優(yōu)化方法的MATLAB實(shí)現(xiàn)〔２學(xué)時〕實(shí)驗(yàn)的目的和要求：通過本次實(shí)驗(yàn)使學(xué)生較為熟練使用MATLAB軟件，并能利用該軟件進(jìn)行約束最優(yōu)化方法的計(jì)算。實(shí)驗(yàn)內(nèi)容：１、罰函數(shù)法的MATLAB實(shí)現(xiàn)２、可行方向法的MATLAB實(shí)現(xiàn)學(xué)習(xí)建議：本次實(shí)驗(yàn)就是要通過對一些具體問題的分析進(jìn)一步熟悉軟件的操作并加深對理論知識的理解。重點(diǎn)和難點(diǎn)：可行點(diǎn)和輔助函數(shù)選取。一罰函數(shù)法用罰函數(shù)法解minf(x)=(x1-1)^2+x2^2S.T.g(x)=x2-1>=0編寫下面m文件fahanshu.msymsx1x2f=(x1-1)^2+x2^2;g=x2-1;x=[x1;x2];x0=[0;0];e=0.0001;M=1;whileabs(subs(g,x,x0))>eifsubs(g,x,x0)<0Q=f+M*g^2;elseQ=f;endx0=xzNewton(Q,x,x0,0.0001);M=10*M;endx0運(yùn)行得：fahanshux0=1.00000.9999與理論值x=[1;1]很接近。二投影梯度法1.梯度投影法根本原理和步驟思想：當(dāng)?shù)c(diǎn)是可行域的內(nèi)點(diǎn)時，將目標(biāo)函數(shù)負(fù)梯度作為搜索方向，當(dāng)?shù)c(diǎn)在可行域邊界上時，將目標(biāo)函數(shù)負(fù)梯度在可行域邊界上的投影作為搜索方向。無論何種情況，所構(gòu)造的方向都是可行下降方向。然后在可行域內(nèi)沿該方向進(jìn)行最優(yōu)一維搜索得到新的迭代點(diǎn)。MATLAB實(shí)現(xiàn):2.代碼及數(shù)值算例:程序源代碼:function[X,FMIN,K]=tidutouying(f,A,b,x1,x,e)%[X,FMIN,K]=tidutouying(f,A,b,x1,e)梯度投影法%f目標(biāo)函數(shù)%A約束矩陣b右端項(xiàng)%x1初始點(diǎn)x自由變量%e精度要求%X極小點(diǎn)%FMIN極小值%K迭代次數(shù)%張超編寫與2014/5/3count=1;n=length(x1);tf=jacobian(f,x)';while1[A1,A2,b1,b2,k]=fenjie(A,b,x1);while1M=A1;ifisempty(M)P=eye(n);elseP=eye(n)-M'*(M*M')^(-1)*M;endPk=-P*subs(tf,x,x1);ifnorm(Pk)<=eifisempty(M)x1;break;elseW=(M*M')^(-1)*M*subs(tf,x,x1);u=W;ifmin(u)<0fori=1:length(u)ifu(i)==min(u)j=i;endendA1(j,:)=[];elsex1;break;endendelseb_=b2-A2*x1;P_=A2*Pk;fori=1:length(P_)ifP_(i)<0r(i)=b_(i)/P_(i);elser(i)=