8.1解優化問題的MATLAB程序結構及其運行機制

8.2線性規劃(LinearProgramming)

8.3二次規劃

8.4無約束非線性規劃(UnconstrainedNonlinearMinimization)

8.5約束非線性規劃

8.6最小二乘規劃

8.7非線性方程求根

第8單元MATLAB優化設計

解優化問題的MATLAB程序具有如圖8-1所示的結構。8.1解優化問題的MATLAB

程序結構及其運行機制圖8-1解優化問題的MATLAB程序結構和數據流主程序至少包括必要數據、初始解、優化設置和解算函數四個部分。必要數據是需預先賦值的優化模型常量參數；初始解是設計變量的初始指定值；優化設置由optimset函數實現，它控制解算函數的算法、設計變量允差、目標函數允差、約束條件允差等；解算函數亦稱優化解法器，是優化工具箱提供的M-file函數，它完成特定的優化模型解法。

優化設置舉例見表8-1。表8-1優化設置語句及其含義線性規劃的MATLAB標準型如下所示，它采用OptimizationToolbox_linprog函數解算。8.2線性規劃(LinearProgramming)

【例8.1】解下面的線性規劃問題

首先將原線性規劃問題轉換為線性規劃的MATLAB標準型，如下所示：其中：

MATLAB解線性規劃程序如下：

clc;clearall;closeall;

C=[-4-5-1];A=[-3-2-1;210];b=[-17;10];Aeq=[111];beq=10;%必要數據

lb=[0;0;0];ub=[]; %必要數據

X0=[0;0;0];

%初始解

options=optimset('TolX',1e-7,'TolFun',1e-7,'TolCon',1e-7); %優化設置

[X,Y,exitflag,output,lambda]=linprog(C,A,b,Aeq,beq,lb,ub,X0,options)%解算程序執行的部分結果如下：

X=

0.0000

10.0000

0.0000

Y=-50.0000

exitflag=1

問題的最優解X和最優值Y分別為

，

【例8.2】某廠生產CP1、CP2、CP3三種產品，均需A、B兩道工序加工。設A工序可在設備A1或A2上完成，B工序可在設備B1、B2或B3上完成。已知產品CP1可在A、B任何一種設備上加工；產品CP2可在任何規格的A設備上加工，但完成B工序時，只能在B1設備上加工；產品CP3只能在A2與B2設備上加工。加工單位產品所需工序時間及其它各項數據見表8-2。試建立問題的線性規劃模型，求解使該廠獲利最大的最優生產規劃。表8-2與產品加工、設備能力和銷售等有關的數據

(1)建立問題的線性規劃模型，步驟如下：

確定設計變量：

x1—A1加工的CP1件數x2—A2加工的CP1件數

x3—B1加工的CP1件數x4—B2加工的CP1件數

x5—B2加工的CP1件數x6—A1加工的CP2件數

x7—A2加工的CP2件數x8—B1加工的CP2件數

x9—A2加工的CP3件數x10—B2加工的CP3件數確定目標函數(總利潤)：確定線性規劃模型：(2)將原線性規劃問題轉換為線性規劃的MATLAB標準型：

(3)?MATLAB解線性規劃程序如下：

clc;clearall;closeall;

C=[-0.75-0.79-0.64-0.56-0.65-1.15-1.38-1.17-1.94-1.09];%必要數據

A=[50000100000;07000090120;0060000800;

00040000011;0000500000]; %必要數據

b=[6000;10000;4000;7000;4000]; %必要數據

Aeq=[11-1-1-100000;0000011-100;000000001-1]; %必要數據

beq=[0;0;0]; %必要數據

lb=[0;0;0;0;0;0;0;0;0;0];ub=[]; %必要數據

X0=[0;0;0;0;0;0;0;0;0;0]; %初始解

options=optimset('TolX',1e-7,'TolFun',1e-7,'TolCon',1e-7); %優化設置

[X,Y,exitflag,output,lambda]=linprog(C,A,b,Aeq,beq,lb,ub,x0,options)%解算

xCP1=X(1)+X(2)

xCP2=X(6)+X(7)

xCP3=X(9)

fval=-Y

(4)程序執行的部分結果如下：

Optimizationterminated.

X=1.0e+003*

1.2000

0.2300

0.0000

0.8586

0.5714

0.0000

0.5000

0.5000

0.3241

0.3241

Y=-4.6017e+003

exitflag=1

xCP1=1.4300e+003產品1總件數

xCP2=500.0000

產品2總件數

xCP3=324.1379

產品3總件數

fval=1.1906e+003總利潤采用OptimizationToolbox_quadprog函數可解算如下所示的MATLAB標準型二次規劃。8.3二次規劃

【例8.3】試求解下面的二次規劃首先將原二次規劃問題轉換為二次規劃的MATLAB標準型，如下所示：

其中：海森矩陣：

采用OptimizationToolbox_quadprog函數編程求解二次規劃。程序如下：

clc;clearall;closeall;

H=[1-1;-12];C=[-2;-6];A=[11;-12;21];b=[2;2;3]; %必要數據

Aeq=[];beq=[];lb=[0;0];ub=[]; %必要數據

X0=[0;0]; %初始解

options=optimset('TolX',1e-7,'TolFun',1e-7,'TolCon',1e-7); %優化設置

[X,Y,exitflag,output,lambda]=quadprog(H,C,A,b,Aeq,beq,lb,ub,X0,options)程序執行的部分結果如下：

Optimizationterminated.

X=最優解

0.6667

1.3333

Y=-8.2222最優值

exitflag=1歐氏空間En、n元自變量X、目標函數f

(X)無約束非線性規劃的MATLAB標準型為：8.4無約束非線性規劃(UnconstrainedNonlinearMinimization)8.4.1無約束非線性規劃的單純形解法

用OptimizationToolbox_fminsearch函數可實現無約束非線性規劃的單純形解法。

【例8.4】求解下面函數的極小點和極小值

采用OptimizationToolbox_fminsearch函數編程求解無約束非線性規劃。程序如下：

clc;clearall;closeall;

fun=@(x)100*(x(2)-x(1)^2)^2+(sqrt(2)-x(1))^2; %目標函數

X0=[0;0]; %初始解

options=optimset('TolX',1e-7,'TolFun',1e-7); %優化設置

[X,Fx,exitflag,output]=fminsearch(fun,X0,options) %解算程序執行的部分結果如下：

X=極小點

1.4142

2.0000

Fx=9.3621e-016極小值

exitflag=18.4.2無約束非線性規劃的梯度解法

用OptimizationToolbox_fminunc函數可實現無約束非線性規劃的梯度解法，它需要調用計算目標函數和梯度的M-file函數，故要先完成目標函數和梯度計算的M-file函數編程。

【例8.5】求解下面的無約束非線性規劃

梯度計算公式：首先為目標函數和梯度計算編寫一個M-file函數uncfun.m。程序如下：

function[f,g]=uncfun(x)

f=3*x(1)^2+2*x(1)*x(2)+x(2)^2;

ifnargout>1%若輸出多于1

g(1)=6*x(1)+2*x(2);

g(2)=2*x(1)+2*x(2);

end

采用OptimizationToolbox_fminunc函數編程調用uncfun并求解無約束非線性規劃。程序如下：X0=[0;0];%初始解

options=optimset('GradObj','on','TolX',1e-7,'TolFun',1e-7); %優化設置

[X1,Fx1,exitflag,output,grad]=fminunc(@uncfun,X0)%缺省設置解算

[X2,Fx2,exitflag,output,grad]=fminunc(@uncfun,X0,options)%用戶設置解算

程序執行的部分結果如下：

X1=1.0e-005*

0.2675

-0.2231

Fx1=1.4509e-011

exitflag=1

grad=1.0e-004*

0.1163

0.0087X2=1.0e-022*

-0.1323

0

Fx2=5.2549e-046

exitflag=1

grad=1.0e-022*

-0.7941

-0.26478.4.3無約束非線性規劃的最速下降解法

用OptimizationToolbox_fminunc函數可實現無約束非線性規劃的最速下降解法。steepdesc最速下降法需通過優化設置啟用。

【例8.6】用最速下降法求解下面的無約束非線性規劃，選定初始解X0?=?[0;0]采用OptimizationToolbox_fminunc函數可實現最速下降法解無約束非線性規劃。程序如下：

clc;clearall;closeall;

X0=[0;0];%初始解

uncfun=@(x)x(1)^2+2*x(2)^2-2*x(1)*x(2)-3*x(2); %目標函數

options=optimset('LargeScale','off','HessUpdate','steepdesc',…

'MaxFunEvals',270); %優化設置

[X,Fx,exitflg]=fminunc(uncfun,X0,options) %用最速下降法解算

程序執行的結果如下：

X=

1.5000

1.5000

Fx=-2.2500

exitflg=18.4.4無約束非線性規劃的DFP變尺度解法

用OptimizationToolbox_fminunc函數可實現DFP擬牛頓變尺度法解無約束非線性規劃。DFP需通過優化設置啟用。

【例8.7】用DFP求解下面的非線性規劃，選定初始解X0?=?[4,0]。采用OptimizationToolbox_fminunc函數可實現DFP擬牛頓法編程求解無約束非線性規劃。程序如下：

clc;clearall;closeall;

X0=[4;0]; %初始解

uncfun=@(x)-1./(x(1)^2+x(2)^2+2); %目標函數

options=optimset('LargeScale','off','HessUpdate','dfp');%優化設置

[X,Fx,exitflg]=fminunc(uncfun,X0,options) %用DFP法解算

程序執行的結果如下：

X=1.0e-006*

-0.7840

0

Fx=-0.5000

exitflg=18.4.5無約束非線性規劃的BFGS變尺度解法

用OptimizationToolbox_fminunc函數可實現BFGS法解無約束非線性規劃。BFGS法需通過優化設置啟用。

【例8.8】用BFGS法求解下面的非線性規劃，選定初始解X0?=?[4,0]。采用OptimizationToolbox_fminunc函數可實現BFGS法編程求解無約束非線性規劃。程序如下：

clc;clearall;closeall;

X0=[4;0]; %初始解

uncfun=@(x)-1./(x(1)^2+x(2)^2+2); %目標函數

options=optimset('LargeScale','off','HessUpdate','bfgs'); %優化設置

[X,Fx,exitflg]=fminunc(uncfun,X0,options) %用BFGS法解算

程序執行的結果如下：

X=1.0e-006*

-0.7840

0

Fx=-0.5000

exitflg=18.5.1約束一元非線性規劃

約束一元非線性規劃求解某區間上一元非線性函數的最小化問題，它的MATLAB標準型為：8.5約束非線性規劃

【例8.9】求解下面的約束一元非線性規劃問題

采用OptimizationToolbox_fminbnd函數編程求解約束一元非線性規劃。程序如下：

clc;clearall;closeall;

objfun=@(x)(x-1.5)^2+0.2*x^3; %目標函數

options=optimset('TolX',1e-7,'TolFun',1e-7,'TolCon',1e-7);%優化設置

[X,Fx,exitflg]=fminbnd(objfun,0,3,options) %解算

程序執行的結果如下：

X=

1.1222

Fx=

0.4254

exitflg=18.5.2約束多元非線性規劃

約束多元非線性規劃問題的MATLAB標準型為：非線性不等式約束非線性等式約束線性不等式約束線性等式約束設計變量區間約束

【例8.10】求解下面的約束多元非線性規劃問題，初始解選定X0?=?[10;10;10]。

首先將約束多元非線性規劃問題轉換為MATLAB標準型。

其中：采用OptimizationToolbox_fmincon函數編程求解約束多元非線性規劃。程序如下：

clc;clearall;closeall;

objfun=@(x)-x(1)*x(2)*x(3); %目標函數X0=[10;10;10]; %初始解A=[-1-2-2;122];b=[0;72]; %必要數據options=optimset(‘TolX’,1e-7,‘TolFun’,1e-7,‘TolCon’,1e-7);%優化設置

[X,Y,exitflg]=fmincon(objfun,X0,A,b,[],[],[],[],[],options)%解算程序執行的結果如下：

Activeinequalities(towithinoptions.TolCon=1e-007):

lowerupperineqlinineqnonlin

2

X=

24.0000

12.0000

12.0000

Y=-3.4560e+003

exitflg=58.5.3約束多元最大最小規劃

最大最小規劃是指實現約束多元非線性向量函數的最大分量最小化，它的MATLAB標準型為：非線性不等式約束非線性等式約束線性不等式約束線性等式約束設計變量區間約束

【例8.11】求解下面的最大最小規劃，初始解選定X0?=?[0.1;0.1]。

其中

例8.11的最大最小化問題無約束，本身已是MATLAB標準型，無需模型轉換。首先將目標函數編寫成M-file函數minmaxobj.m。程序如下：

functionf=minmaxobj(x);

f(1)=2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304;

f(2)=-x(1)^2-3*x(2)^2;

f(3)=x(1)+3*x(2)-18;

f(4)=-x(1)-x(2);

f(5)=x(1)+x(2)-8;采用OptimizationToolbox_fminimax函數編程調用minmaxobj并求解最大最小規劃。程序如下：

clc;clearall;closeall;

x0=[0.1;0.1]; %初始解

options=optimset('TolX',1e-7,'TolFun',1e-7,'TolCon',1e-7); %優化設置

formatshorte;

[x,fval,maxfval,exitflag]=fminimax(@minmaxobj,x0,[],[],[],[],[],[],[],options)

程序執行的結果如下：

Activeinequalities(towithinoptions.TolCon=1e-007):%起作用約束

lowerupperineqlinineqnonlin

1

5

x= %最優解

4.0000e+000

4.0000e+000

fval=

4.9454e-011-6.4000e+001-2.0000e+000-8.0000e+000-1.5445e-012 %目標函數值

maxfval=4.9454e-011 %最大目標函數值

exitflag=48.5.4約束多元多目標規劃

多目標規劃問題是讓非線性向量函數(多目標)整體上與期望目標的差異最小，它的MATLAB標準型為：非線性不等式約束非線性等式約束線性不等式約束線性等式約束設計變量區間約束

【例8.12】例8.11最大最小規劃問題采用多目標規劃模型求解，初始解仍選定X0?=?[0.1;0.1]。

確定目標函數向量F(X)的目標值向量goal和目標權重向量weight：

；先將目標函數編寫成M-file函數multgoalobj.m。程序如下：

functionf=multgoalobj(x);

f(1)=2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304;

f(2)=-x(1)^2-3*x(2)^2;

f(3)=x(1)+3*x(2)-18;

f(4)=-x(1)-x(2);

f(5)=x(1)+x(2)-8;采用OptimizationToolbox_fgoalattain函數編程調用multgoalobj并求解多目標規劃。程序如下：

clc;clearall;closeall;

x0=[0.1;0.1]; %初始解

goal=[00000]'; %目標值

weight=goal+1; %目標權重，選定權重相同

options=optimset('TolX',1e-7,'TolFun',1e-7,'TolCon',1e-7);%優化設置

formatshorte;

[x,fval,attainfactor,exitflag]=fgoalattain(@minmaxobj,x0,goal,weight,...

[],[],[],[],[],[],[],options)

程序執行的結果如下：

Activeinequalities(towithinoptions.TolCon=1e-007):%起作用約束

lowerupperineqlinineqnonlin

1

5

x=%最優解

4.0000e+000

4.0000e+000

fval=5.2978e-011-6.4000e+001-2.0000e+000-8.0000e+000-1.6502e-012

attainfactor=-1.7365e-012%松弛因子值

exitflag=58.6.1非線性方程最小二乘規劃

響應變量的測定值Ydata與它的估計值f?(X)之差稱做誤差，自變量X的n點測試就會產生n個誤差。采用OptimizationToolbox_lsqcurvefit函數編程實現非線性方程最小二乘規劃，即實現n個誤差平方和的最小化，它的MATLAB標準型為：8.6最小二乘規劃其中

響應變量Y與自變量X的關系可表示為下面的響應模型：維模型參數，包括β0；維自變量樣本(試驗數據)；維響應變量樣本(試驗數據)；維響應變量估計值。是1維響應變量是p維自變量是p+1維模型參數是誤差

【例8.13】為考察水稻品種IR72的籽粒平均粒重y?(毫克)與開花后天數x的關系，測定了兩變量的9對數據，詳見表8-3。表8-3水稻品種IR72的性狀觀測為試驗測定數據擬合一個Logistic生長模型，即：

，(參數k、a、b均大于零)

采用lsqcurvefit函數的其他案例，可考察第4單元中4.4.2小節的一元函數最小二乘擬合問題和第7單元中7.7.1小節的一元非線性回歸問題。樣本數據和模型參數表為下面的向量：Logistic生長模型對參數β的偏導數表示為梯度G的轉置向量：

n個試驗觀測上的偏導數向量表示為Jacobian矩陣：將目標函數和Jacobian矩陣編寫成M-file函數IR72objfun.m。程序如下：

function[F,J]=IR72objfun(beta,xdata)

fenmu=1+beta(2)*exp(-beta(3)*xdata);

F=beta(1)./fenmu;

ifnargout>1;%Jacobianofthefunctionevaluatedatbeta

yk=1./fenmu;

ya=-beta(1)*exp(-beta(3)*xdata)./fenmu.^2;

yb=beta(1)*beta(2)*xdata.*exp(-beta(3)*xdata)./fenmu.^2;

J=[ykyayb];

end;采用OptimizationToolbox_lsqcurvefit函數編程調用IR72objfun并求解最小二乘規劃。程序如下：

clc;clearall;

xdata=[03691215182124]';

ydata=[0.30.723.319.7113.0916.8517.7918.2318.43]';

beta0=[1,1,1]';%Startingguess

lb=[000]';ub=[];

options=optimset('GradObj','on','TolX',1e-7,'TolFun',1e-7,'TolCon',1e-7);

[beta,resnorm,residual,exitflag]=...

lsqcurvefit(@IR72objfun,beta0,xdata,ydata,lb,ub,options)程序執行的部分結果如下：

Optimizationterminated:relativefunctionvaluechangingbylessthanOPTIONS.TolFun.

beta=

1.8280e+001

4.8221e+001

4.1987e-001

resnorm=2.2398e+000

exitflag=3

平均粒重與開花后天數關系的Logistic生長模型擬合如圖8-2所示。圖8-2平均粒重與開花后天數關系的Logistic生長模型擬合8.6.2非線性方程組最小二乘規劃

采用OptimizationToolbox_lsqnonlin函數可編程實現非線性方程組或向量函數的最小二乘規劃，問題的MATLAB標準型為：

其中，為非線性函數組或向量函數。

【例8.14】求解下面的非線性方程組最小二乘規劃：

其中Jacobian矩陣的計算：將目標函數和Jacobian矩陣的計算編寫成M-file函數nonlinobjfun.m。程序如下：

function[FJ]=nonlinobjfun(x);

i=[1:10]’;

F=2+2*i-exp(i*x(1))-exp(i*x(2));

ifnargout>1;%twooutputarguments

fori=1:10;

forj=1:2;

J(i,j)=-i*exp(i*x(j));

end;

end;%Jacobianofthefunctionevaluatedatx

end;采用OptimizationToolbox_lsqnonlin函數編程調用nonlinobjfun并求解最小二乘規劃。程序如下：

clc;clearall;

options=optimset('Jacobian','on','TolX',1e-7,'TolFun',1e-7);

x0=[0.3;0.4];%Startingguess

[x,resnorm,residual,exitflag]=lsqnonlin(@nonlinobjfun,x0,[],[],options)

程序執行的部分結果如下：

Optimizationterminated:normofthecurrentstepislessthanOPTIONS.TolX.

x=

0.2578

0.2578

resnorm=

124.3622

exitflag=28.6.3線性方程組約束最小二乘規劃

采用OptimizationToolbox_lsqlin函數可編程實現線性方程組約束最小二乘規劃，常用于求取約束或無約束超定線性方程組的最小二乘解。問題的MATLAB標準型為：

【例8.15】求解下面的線性方程組約束最小二乘規劃：

其中采用OptimizationToolbox_lsqlin函數編程求解線性方程組約束最小二乘規劃。程序如下：

clc;clearall;

C=[0.95010.76200.61530.4057

0.23110.45640.79190.9354

0.60680.01850.92180.9169

0.48590.82140.73820.4102

0.89120.44470.17620.8936];

d=[0.05780.35280.81310.00980.1388]’;

lb=-0.1*ones(4,1);ub=2*ones(4,1);

A=[0.20270.27210.74670.4659

0.19870.19880.44500.4186

0.60370.01520.93180.8462];

b=[0.52510.20260.6721]’;

options=optimset('LargeScale','off','TolX',1e-7,'TolFun',1e-7,'TolCon',1e-7);

[x,resnorm,residual,exitflag]=lsqlin(C,d,A,b,[],[],lb,ub,[],options)

程序執行的部分結果如下：

Optimizationterminated.

x=

-0.1000

-0.1000

?0.2152

?0.3502

resnorm=0.1672

exitflag=18.6.4線性方程組非負最小二乘規劃

采用OptimizationToolbox_lsqnonneg函數可編程實現線性方程組非負最小二乘規劃，常用于求取超定線性方程組的非負約束最小二乘解。問題的MATLAB標準型為：

【例8.16】求下面二元線性方程組的非負最小二乘解，即求的最小二乘解。

其中采用OptimizationToolbox_lsqnonneg函數編程求解線性方程組非負最小二乘規劃。程序如下：

clc;clearall;

C=[0.03720.2869

0.68610.7071

0.62330.6245

0.6344