非線性最小二乘數(shù)據(jù)擬合高斯牛頓法

非線性最小二乘數(shù)據(jù)擬合高斯牛頓法第1頁(yè)，共33頁(yè)，2023年，2月20日，星期四例題1:某篩分試驗(yàn)測(cè)得篩孔尺寸與累計(jì)產(chǎn)率間的關(guān)系如表所示，試建立模型。

產(chǎn)率與篩孔關(guān)系篩孔/mm502513630.5產(chǎn)率/%1.51130557088此為一近似指數(shù)關(guān)系，即型。

第2頁(yè)，共33頁(yè)，2023年，2月20日，星期四用MATLAB編程首先建立一個(gè)′ep′函數(shù)文件。functionf=ep(b,x)；f=(b(1)*exp(－b(2)*x)第3頁(yè)，共33頁(yè)，2023年，2月20日，星期四輸入x=[502513630.5]’;y=[1.51130557088];b0=[11]; %初始值[b,r,j]=nlinfit(x,y,′ep′,b0);b[yp,d]=nlperdci(′ep′,x,b,r,j);ci=nlparci(b,r,j)%輸出預(yù)測(cè)值、殘差及置信區(qū)間nlintool(x,y,′ep′,b0); %繪制非線性擬合曲線圖3.14plot(x,r,′――′) %繪殘差曲第4頁(yè)，共33頁(yè)，2023年，2月20日，星期四結(jié)果得b=91.41580.0855ci=90.208592.62310.08280.0883移動(dòng)垂直于x軸的虛線可以顯示y隨x的變化。數(shù)學(xué)模型為y=91.4158e-0.0855x第5頁(yè)，共33頁(yè)，2023年，2月20日，星期四例2已知浮煤累計(jì)灰分及產(chǎn)率為x和y，試進(jìn)行非線性擬合。x=[3.467.339.4710.4811.1911.7912.7820.50]；y=[10.6956.8476.8982.1584.7086.3288.45100]；現(xiàn)采用含有雙曲正切的函數(shù)擬合。f=100*(b(1)+b(2)*x+b(2)*b(3)*tanh(b(4)*(x-b(5))))；

第6頁(yè)，共33頁(yè)，2023年，2月20日，星期四輸入x=x’；b0=[0.20-0.4-0.15]；%初始值[b,r,j]=nlinfit(x,y,′mf′,b0)；[yp,d]=nlpredci(′mf′,x,b,r,j);nlintool(x,y,′mf′,b0)第7頁(yè)，共33頁(yè)，2023年，2月20日，星期四結(jié)果b=0.34140.0136-0.3797-0.33656.2948即y=100*(0.3414+0.0136x-0.3797tanh(-0.3365(x-6.2948)))由計(jì)算出的yy可知，擬合效果很好。第8頁(yè)，共33頁(yè)，2023年，2月20日，星期四第三講線性規(guī)劃

線性規(guī)劃問題是目標(biāo)函數(shù)和約束條件均為線性函數(shù)的問題minsub.to：

其中f、x、b、beq、lb、ub為向量，A、Aeq為矩陣。其它形式的線性規(guī)劃問題都可經(jīng)過適當(dāng)變換化為此標(biāo)準(zhǔn)形式。第9頁(yè)，共33頁(yè)，2023年，2月20日，星期四補(bǔ)充知識(shí):解方程或方程組用solve符號(hào)解代數(shù)方程solve(′eq1′，′eq2′,′eqn′,′var1,var2,…,varn′)eqn為符與方程，varn為符號(hào)變量。

解代數(shù)方程ax2+bx+c=0的根。symsabcxg=a*x^2+b*x+c；s=solve(g)s=[1/2/a*(-g+(g^2-4*2*c)^(1/2))][1/2/a*(-g+(g^2+4*2*c)^(1/2))]類似地，solve(′2*x^2+6*x+4′)得ans=[-2][-1]第10頁(yè)，共33頁(yè)，2023年，2月20日，星期四求解x2+xy+y=3，x為參數(shù)。solve(′x^2+x*y+y=3′,′x′)ans=[-1/2*y+1/2*(y^2-4*y+12)^(1/2)][-1/2*y-1/2*(y^2-4*y+12)^(1/2)]第11頁(yè)，共33頁(yè)，2023年，2月20日，星期四解au2+υ2=0,u-υ=1。[u,v]=solve(′a*u^2+v^2=0′,′u-v=1′)u=[1/2/(a+1)*(-2*a+2*(-a)^(1/2))+1][1/2/(a+1)*(-2*a-2*(-a)^(1/2))+1]v=[1/2/(a+1)8(-2*a+2*(-a)^(1/2))][1/2/(a+1)8(-2*a-2*(-a)^(1/2))]第12頁(yè)，共33頁(yè)，2023年，2月20日，星期四2解線性方程組MATLAB用函數(shù)linsolve求解線性方程組Ax=b，要求A的列數(shù)等于b的行數(shù)；也可用矩陣除等方法求解。linsolve的語(yǔ)法格式為x=linsolve(A,b)第13頁(yè)，共33頁(yè)，2023年，2月20日，星期四例題解線性方程組第14頁(yè)，共33頁(yè)，2023年，2月20日，星期四a=[121；223；-1-30]；%rank(a)=3，有惟一解b=[032]’;x=linsolve(a,b)x=1-11第15頁(yè)，共33頁(yè)，2023年，2月20日，星期四3非線性方程

非線性方程的標(biāo)準(zhǔn)形式為f(x)=0函數(shù)

fzero格式x=fzero(fun,x0)%用fun定義表達(dá)式f(x)，x0為初始解。x=fzero(fun,x0,options)[x,fval]=fzero(…)%fval=f(x)[x,fval,exitflag]=fzero(…)[x,fval,exitflag,output]=fzero(…)說明該函數(shù)采用數(shù)值解求方程f(x)=0的根。第16頁(yè)，共33頁(yè)，2023年，2月20日，星期四非線性方程解例題求x^3-2*x-5=0的根解：>>fun='x^3-2*x-5';>>z=fzero(fun,2)%初始估計(jì)值為2結(jié)果為z=2.0946第17頁(yè)，共33頁(yè)，2023年，2月20日，星期四非線性方程組求解非線性方程組的標(biāo)準(zhǔn)形式為：F(x)=0其中：x為向量，F(xiàn)(x)為函數(shù)向量。函數(shù)

fsolve格式x=fsolve(fun,x0)%用fun定義向量函數(shù)，其定義方式為：先定義方程函數(shù)functionF=myfun(x)。F=[表達(dá)式1；表達(dá)式2；…表達(dá)式m]%保存為myfun.m，并用下面方式調(diào)用：x=fsolve(@myfun,x0)，x0為初始估計(jì)值。x=fsolve(fun,x0,options)[x,fval]=fsolve(…)%fval=F(x)，即函數(shù)值向量[x,fval,exitflag]=fsolve(…)[x,fval,exitflag,output]=fsolve(…)[x,fval,exitflag,output,jacobian]=fsolve(…)%jacobian為解x處的Jacobian陣。其余參數(shù)與前面參數(shù)相似。第18頁(yè)，共33頁(yè)，2023年，2月20日，星期四非線性方程組求解例題解下列非線性方程組第19頁(yè)，共33頁(yè)，2023年，2月20日，星期四程序及做法:functiony=xzz(x)；y(1)=x(1)+2*x(2)+x(3);y(2)=2*x(1)+2*x(2)+3*x(3)-3;y(3)=-x(1)-3*x(2)-2;y=[y(1)y(2)y(3)];x0=[111];[x,fva1,exitflag,output]=fsolve(′xzz′,x0)第20頁(yè)，共33頁(yè)，2023年，2月20日，星期四x=1.0000-1.00001.0000fval=1.0e-006*-0.27290.0810-0.1294exitflag=1output=firstorderopt:2.9957e-008iterations:4func-count:17cgiterations:3algorithm:[1x43char]第21頁(yè)，共33頁(yè)，2023年，2月20日，星期四求解下列非線性方程組x(1)-0.7*sin(x(1))-0.2*cos(x(2))=0x(2)-0.7*sin(x(1))+0.2*cos(x(2))=0函數(shù)文件functionf=xzz(x);f(1)=x(1)-0.7*sin(x(1))-0.2*cos(x(2));f(2)=x(2)-0.7*sin(x(1))+0.2*cos(x(2));f=[f(1)f(2)];輸入x0=[11]；x=fsolve(′xzz′，x0)x=0.52650.5079第22頁(yè)，共33頁(yè)，2023年，2月20日，星期四線性規(guī)劃問題是目標(biāo)函數(shù)和約束條件均為線性函數(shù)的問題minsub.to：

其中f、x、b、beq、lb、ub為向量，A、Aeq為矩陣。其它形式的線性規(guī)劃問題都可經(jīng)過適當(dāng)變換化為此標(biāo)準(zhǔn)形式。第23頁(yè)，共33頁(yè)，2023年，2月20日，星期四函數(shù)linprog格式x=linprog(f,A,b)%求minf'*xsub.to線性規(guī)劃的最優(yōu)解。x=linprog(f,A,b,Aeq,beq)%等式約束，若沒有不等式約束，則A=[]，b=[]。x=linprog(f,A,b,Aeq,beq,lb,ub)%指定x的范圍，若沒有等式約束，則Aeq=[]，beq=[]x=linprog(f,A,b,Aeq,beq,lb,ub,x0)%設(shè)置初值x0x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options)%options為指定的優(yōu)化參數(shù)[x,fval]=linprog(…)%返回目標(biāo)函數(shù)最優(yōu)值，即fval=f'*x。[x,lambda,exitflag]=linprog(…)%lambda為解x的Lagrange乘子。[x,lambda,fval,exitflag]=linprog(…)%exitflag為終止迭代的錯(cuò)誤條件。[x,fval,lambda,exitflag,output]=linprog(…)%output為關(guān)于優(yōu)化的一些信息說明若exitflag>0表示函數(shù)收斂于解x，exitflag=0表示超過函數(shù)估值或迭代的最大數(shù)字，exitflag<0表示函數(shù)不收斂于解x；若lambda=lower表示下界lb，lambda=upper表示上界ub，lambda=ineqlin表示不等式約束，lambda=eqlin表示等式約束，lambda中的非0元素表示對(duì)應(yīng)的約束是有效約束；output=iterations表示迭代次數(shù)，output=algorithm表示使用的運(yùn)算規(guī)則，output=cgiterations表示PCG迭代次數(shù)。第24頁(yè)，共33頁(yè)，2023年，2月20日，星期四求下面的優(yōu)化問題minsub.to第25頁(yè)，共33頁(yè)，2023年，2月20日，星期四f=[-5;-4;-6];A=[1-11;324;320];b=[20;42;30];lb=zeros(3,1);[x,fval,exitflag,output,lambda]=linprog(f,A,b,[],[],lb)第26頁(yè)，共33頁(yè)，2023年，2月20日，星期四結(jié)果為：x=%最優(yōu)解

0.000015.00003.0000fval=%最優(yōu)值

-78.0000exitflag=%收斂

1output=iterations:6%迭代次數(shù)

cgiterations:0algorithm:'lipsol'%所使用規(guī)則lambda=ineqlin:[3x1double]eqlin:[0x1double]upper:[3x1double]lower:[3x1double]>>lambda.ineqlinans=0.00001.50000.5000>>lambda.lowerans=1.00000.00000.0000第27頁(yè)，共33頁(yè)，2023年，2月20日，星期四某廠有三個(gè)礦的來煤。各礦煤的煤質(zhì)特征及要求的配煤指標(biāo)列入表5.1中。試確定怎樣摻配才能達(dá)到配煤各項(xiàng)指標(biāo)的要求，同時(shí)又能使原料煤成本最低（不含加工費(fèi)）。表5.1三種單煤及要求的配煤指標(biāo)指標(biāo)Mt/%Ad/%Vdaf/%CRCst,d/%Qnet,ar/kcal·kgHGIST/C成本/元·t-1甲煤7.7820.4317.0440.335974851400255.00乙煤10.8725.4133.6242.325177541290244.30丙煤3.0525.6334.0430.475515501450220.75要求指標(biāo)≤7.00≥23.0≤5≤0.80≥5500≥55≥1300最低第28頁(yè)，共33頁(yè)，2023年，2月20日，星期四設(shè)甲、乙、丙三種煤的配比分別為x1,x2,x3。根據(jù)指標(biāo)的可加性，利用線性規(guī)劃來計(jì)算各最佳配煤參數(shù)。列出目標(biāo)函數(shù)及各指標(biāo)的約束方程。目標(biāo)函數(shù)minf(x)=255x1+244.3x2+220.75x3約束條件7.78x1+10.87x2+3.05≤7.00%水分17.04x1+33.62x2+34.04x3≥23.00 %揮發(fā)分4x1+4x2+3x3≤5 %焦渣特征0.33x1+2.32x2+0.47x3≤0.8 %全硫5974x1+5177x2+5515x3≥5500 %發(fā)熱量85x1+54x2+50x3≥55 %可磨性指數(shù)1400x1+1290x2+1450x3≥1300 %灰熔融性x1+x2+x3=10≤x1,0≤x2,x3≤0第29頁(yè)，共33頁(yè)，2023年，2月20日，星期四按MATLAB語(yǔ)法輸入各系數(shù)值。f=[255;244.3;220.75];A=[7,7810,873.05;-17.04-33.62-34.04;443;0.332.320.47;-5974-5177-5515;-85-54-50;-1400-1290-1450];b=[7;-23;5;0.80;-5500;-55;-1300];a1=[111];b1=1;lb=zeros(3,1);[x,fval]=linprog(f,A,b,a1,b1,lb)計(jì)算得，x1=0.1429,x2=0,x3=0.8571，成本為225.64元/t。即三種煤的配比為14.29%，0%和85.71%?？紤]到皮帶秤的精度可取整數(shù)。為充分利用各煤種，也可以令三種煤都摻配，這時(shí)需要適當(dāng)調(diào)整各配比（可用lb控制）。本例中三種煤的配比可以是13%，10%，77%，成本227.61元/t，或者13%，11%，76%，成本227.79元/t。但若乙丙或