管理數(shù)學(xué)實(shí)驗(yàn)- MATLAB在管理運(yùn)籌中的應(yīng)用 3-管理數(shù)學(xué)實(shí)驗(yàn)-規(guī)劃論及MATLAB計(jì)算

第3章規(guī)劃論及MATLAB計(jì)算1、線性規(guī)劃2、MATLAB優(yōu)化工具箱簡介3、非線性最優(yōu)化4、動(dòng)態(tài)規(guī)劃5、GUI優(yōu)化工具1、線性規(guī)劃1.1線性規(guī)劃模型1.2用MATLAB優(yōu)化工具箱解線性規(guī)劃1.3應(yīng)用案例1.1線性規(guī)劃的模型（1）問題的提出例1.某工廠在計(jì)劃期內(nèi)要安排甲、乙兩種產(chǎn)品的生產(chǎn)，已知生產(chǎn)單位產(chǎn)品所需的資源A、B、C的消耗以及資源的計(jì)劃期供給量，如下表：問題：工廠應(yīng)分別生產(chǎn)多少單位甲、乙產(chǎn)品才能使工廠獲利最多？項(xiàng)目甲產(chǎn)品乙產(chǎn)品資源供給量資源A23180資源B32210資源C15250單位產(chǎn)品獲利6060

線性規(guī)劃的模型解：設(shè)甲、乙產(chǎn)品的產(chǎn)量分別為x1、x2

，工廠獲利為z，則目標(biāo)函數(shù)：

Maxz=60x1

+60x2

約束條件：s.t.

2x1

+3x2≤180

3x1+2x2≤210

x1+5x2≤250x1,x2≥0甲產(chǎn)品x1乙產(chǎn)品x2資源供給量資源A23180資源B32210資源C15250單位產(chǎn)品獲利6060

線性規(guī)劃的模型從上述的例子看出，建立數(shù)學(xué)模型的基本過程是：1）搞清要解決的問題：目標(biāo)和條件;2）設(shè)置決策變量--描述解決問題的方案;3）描述約束條件和非負(fù)約束;4）給出目標(biāo)函數(shù)，確定目標(biāo)函數(shù)的優(yōu)化方向，即優(yōu)化是對(duì)目標(biāo)函數(shù)取最大還是最小。線性規(guī)劃的模型（2）線性規(guī)劃模型：一般形式目標(biāo)函數(shù)：

Max（Min）z=c1x1+c2x2+…+cnxn約束條件：s.t.

a11x1+a12x2+…+a1nxn≤（=,≥）b1

a21x1+a22x2+…+a2nxn≤（=,≥）b2…………

am1x1+am2x2+…+amnxn≤（=,≥）bm

x1，x2，…，xn≥0線性規(guī)劃的模型線性規(guī)劃一般數(shù)學(xué)模型的矩陣形式目標(biāo)函數(shù)max(或min)z=c·x約束條件A·x≤（=,≥）

bx≥0式中c=(c1,c2,…,cn),x=(x1,x2,…,xn)τa11a12

…a1nA=a21a22

…a2n,b=(b1,b2,…,bm)τ

…

am1am2

…amn

線性規(guī)劃的模型（3）線性規(guī)劃模型：標(biāo)準(zhǔn)形式目標(biāo)函數(shù)：

Maxz=c1x1+c2x2+…+cnxn約束條件：s.t.

a11x1+a12x2+…+a1nxn=b1

a21x1+a22x2+…+a2nxn=b2…………

am1x1+am2x2+…+amnxn=bm

x1，x2，…，xn≥0特點(diǎn)：目標(biāo)最大化約束為等式右端項(xiàng)非負(fù)決策變量非負(fù)線性規(guī)劃的模型線性規(guī)劃標(biāo)準(zhǔn)型的矩陣形式目標(biāo)函數(shù)maxz=c·x約束條件A·x=bx≥0式中c=(c1,c2,…,cn),x=(x1,x2,…,xn)τa11a12

…a1nA=a21a22

…a2n,b=(b1,b2,…,bm)τ

…

am1am2

…amn

一般模型化為標(biāo)準(zhǔn)型的處理過程（4）線性規(guī)劃的圖解法目標(biāo)函數(shù)：

maxZ=X1+X2約束條件：最優(yōu)解：X1=50,X2=6.67可行域可行解最優(yōu)解最優(yōu)值Z=56.67Z=10Z=48線性規(guī)劃的模型線性規(guī)劃的圖解法例1的圖解目標(biāo)函數(shù)：

Maxz=60x1+60x2約束條件：2x1+3x2≤180(A)

3x1+2x2≤210(B)

x1+5x2≤250(C)x1，

x2≥0

(D)得到最優(yōu)解：

x1=54，x2=24

最優(yōu)目標(biāo)值z(mì)=4680解的幾種情況：線性規(guī)劃的最優(yōu)解如果存在則必定有一個(gè)頂點(diǎn)（極點(diǎn)）是最優(yōu)解①唯一解目標(biāo)函數(shù)等值線與約束邊界只有一個(gè)交點(diǎn)②無窮多最優(yōu)解目標(biāo)函數(shù)等值線與約束邊界平行③無界解可行域不封閉④無可行解可行域?yàn)榭占?）線性規(guī)劃解的概念引入松馳變量____含義是資源的剩余量例1中引入s1，s2，s3模型化為標(biāo)準(zhǔn)型目標(biāo)函數(shù)：Maxz=60x1+60x2+0s1+0s2+0s3約束條件：s.t.2x1+3x2+s1=180

3x1+2x2+s2=210

x1+5x2+s3

=

250x1,x2,s1,s2,s3≥0對(duì)于標(biāo)準(zhǔn)型的最優(yōu)解

x1=54x2=24，s1=0s2=0s3=76說明：生產(chǎn)54單位甲產(chǎn)品和24單位乙產(chǎn)品將消耗完所有的A、B資源，但對(duì)資源C則還剩余76。線性規(guī)劃問題解的概念基最優(yōu)解、最優(yōu)解、基可行解、基解、可行解的關(guān)系如下所示：基最優(yōu)解基可行解可行解最優(yōu)解基解線性規(guī)劃的幾何意義和基本定理（6）線性規(guī)劃的基本定理①線性規(guī)劃問題的所有可行解構(gòu)成的集合(可行域)R={x|A·x=b,x≥0}R是一凸集(包括無界域)，它有有限個(gè)頂點(diǎn)；②線性規(guī)劃問題的每個(gè)基可行解對(duì)應(yīng)可行域凸集R的一個(gè)頂點(diǎn)；③若線性規(guī)劃問題有最優(yōu)解,則必定在某頂點(diǎn)處得到基本定理把可行域的有限個(gè)頂點(diǎn)與基可行解這一代數(shù)概念聯(lián)系起來，可通過求基可行解的代數(shù)方法來得到可行域的一切極點(diǎn)，能在有限的計(jì)算中獲得最優(yōu)點(diǎn)。1.2用MATLAB優(yōu)化工具箱解線性規(guī)劃（1）線性規(guī)劃模型1函數(shù)：minz=cx命令：x=linprog（c，A，b）（2）線性規(guī)劃模型2函數(shù)：minz=cx

命令：x=linprog（c，A，b，Aeq,beq）注意：若沒有不等式：存在，則令A(yù)=[]，b=[].（3）線性規(guī)劃模型3函數(shù)：minz=CxVLB≤X≤VUB命令：[1]x=linprog（c，A，b，Aeq,beq,VLB，VUB）

[2]x=linprog（c，A，b，Aeq,beq,VLB，VUB,X0）

注意：[1]若沒有等式約束:,

則令A(yù)eq=[],beq=[].[2]其中x0表示初始點(diǎn)[x,fval]=linprog(…)[x,fval,exitflag]=linprog(…)[x,fval,exitflag,output]=linprog(…)[x,fval,exitflag,output,lambda]=linprog(…)求解線性規(guī)劃模型的輸出擴(kuò)展函數(shù)調(diào)用形式如下：（4）函數(shù)（命令）：[x,fval]=linprog(…)返回最優(yōu)解ｘ及ｘ處的目標(biāo)函數(shù)值fval.x=bintprog(c)x=bintprog(c,A,b)x=bintprog(c,A,b,Ae,be)x=bintprog(c,A,b,Ae,be,x0)x=bintprog(c,A,b,Ae,be,x0,options)[x,fval]=bintprog(…)[x,fval,exitflag]=bintprog(…)[x,fval,exitflag,output]=bintprog(…)（5）0-1整數(shù)規(guī)劃：

minz=cx

或

或求解0-1整數(shù)規(guī)劃模型的函數(shù)調(diào)用形式如下：x=intlinprog(c,intcon,A,b)x=intlinprog(c,intcon,A,b,Ae,be)x=intlinprog(c,intcon,A,b,Ae,be,vlb,vub)x=intlinprog(c,intcon,A,b,Ae,be,vlb,vub,x0)x=intlinprog(c,intcon,A,b,Ae,be,vlb,vub,x0,options)[x,fval]=intlinprog(…)[x,fval,exitflag]=intlinprog(…)[x,fval,exitflag,output]=intlinprog(…)（6）整數(shù)規(guī)劃：

minz=cx

或

或求解整數(shù)規(guī)劃模型的函數(shù)調(diào)用形式如下：解編寫M文件如下：c=[-0.4-0.28-0.32-0.72-0.64-0.6];A=[0.010.010.010.030.030.03;0.02000.0500;00.02000.050;000.03000.08];b=[850;700;100;900];Aeq=[];beq=[];vlb=[0;0;0;0;0;0];vub=[];[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)解:編寫M文件如下：

c=[634];A=[010];b=[50];Aeq=[111];beq=[120];vlb=[30,0,20];vub=[];[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)例3

：某木材公司經(jīng)營的木材貯存在倉庫中，最大貯存量為30萬m3。由于木材價(jià)格隨季節(jié)變化，該公司于每季初購進(jìn)木材，一部分當(dāng)季出售，一部分貯存以后出售。貯存費(fèi)為：a+bu，其中a=90元／m3，b=100元／m3／季，u為貯存的季度數(shù)。由于木材久貯易損，因此當(dāng)年所有庫存木材應(yīng)于秋末售完。各季木材單價(jià)及銷量如下表所示。為獲全年最大利潤，該公司各季應(yīng)分別購銷多少木材?

1.3線性規(guī)劃建模案例各季木材單價(jià)及銷量季節(jié)購進(jìn)價(jià)（元/m3）售出價(jià)（元/m3）最大銷售量（萬m3）冬6200642015春6500666021夏6960704030秋6800688024季節(jié)冬春夏秋max貯存量冬220460-190=270840-290=550680-390=29030春-160540-190=350380-290=9030夏--80-80-190=-27030秋---8030max銷量15213024

假設(shè)：第i季度購買第j季度銷售的木材量為xij（萬m3），則第i季度購買第j季度銷售的木材量為xij（萬m3）的利潤（元/m3）見下表。

各季木材單價(jià)及銷量問題求解的數(shù)學(xué)模型如下：

編寫求解的m文件如下c=-[2202705502901603509080-27080];A=[0111000000;0011011000;0001001010];b=[30;30;30];Aeq=[1000000000;0100100000;0010010100;0001001011];beq=[15;21;30;24];vlb=[0;0;0;0;0;0;0;0;0;0];vub=[];[x,fval]=linprog(-c,A,b,Aeq,beq,vlb,vub)x=[15,0,30,0,21,0,0,0,0,24]’fval=25080優(yōu)化結(jié)果如下：2Matlab優(yōu)化工具箱簡介（1）MATLAB求解優(yōu)化問題的主要函數(shù)（2）優(yōu)化函數(shù)的輸入變量

使用優(yōu)化函數(shù)或優(yōu)化工具箱中其它優(yōu)化函數(shù)時(shí),輸入變量見下表:（3）優(yōu)化函數(shù)的輸出變量（4）控制參數(shù)options的設(shè)置③MaxIter:允許進(jìn)行迭代的最大次數(shù),取值為正整數(shù).Options中常用的幾個(gè)參數(shù)的名稱、含義、取值如下:①Display:顯示水平。取值為'off'時(shí),不顯示輸出;取值為'iter'時(shí),顯示每次迭代的信息;取值為'final'時(shí),顯示最終結(jié)果.默認(rèn)值為'final'。②MaxFunEvals:允許進(jìn)行函數(shù)評(píng)價(jià)的最大次數(shù),取值為正整數(shù).例：opts=optimset('Display','iter','TolFun',1e-8)

該語句創(chuàng)建一個(gè)稱為opts的優(yōu)化選項(xiàng)結(jié)構(gòu),其中顯示參數(shù)設(shè)為

'iter',TolFun參數(shù)設(shè)為1e-8。參數(shù)options可以通過函數(shù)optimset創(chuàng)建或修改。①options=optimset('optimfun')

創(chuàng)建一個(gè)含有所有參數(shù)名,并與優(yōu)化函數(shù)optimfun相關(guān)的默認(rèn)值的選項(xiàng)結(jié)構(gòu)options.②options=optimset('param1',value1,'param2',value2,...)

創(chuàng)建一個(gè)名稱為options的優(yōu)化選項(xiàng)參數(shù),其中指定的參數(shù)具有指定值,所有未指定的參數(shù)取默認(rèn)值.③options=optimset(oldops,'param1',value1,'param2',value2,...)

創(chuàng)建名稱為oldops的參數(shù)的拷貝,用指定的參數(shù)值修改oldops中相應(yīng)的參數(shù).命令的格式如下：3、非線性最優(yōu)化問題3.1無約束優(yōu)化問題3.2非線性規(guī)劃3.1無約束最優(yōu)化問題3.1.1一元函數(shù)無約束優(yōu)化問題3.1.2多元函數(shù)無約束優(yōu)化問題3.1.1一元函數(shù)無約束優(yōu)化問題其中③、④、⑤的等式右邊可選用①或②的等式右邊。函數(shù)fminbnd的算法基于黃金分割法和二次插值法，它要求目標(biāo)函數(shù)必須是連續(xù)函數(shù)，并可能只給出局部最優(yōu)解。常用格式如下：①x=fminbnd(fun,x1,x2)②x=fminbnd(fun,x1,x2

，options)③[x，fval]=fminbnd（...）④[x，fval，exitflag]=fminbnd（...）⑤[x，fval，exitflag，output]=fminbnd（...）主程序?yàn)閙ain1.m:fun=inline('2*exp(-x).*sin(x)');fplot(fun,[0,8]);holdon[xmin,ymin]=fminbnd(fun,0,8);fun1=inline('-2*exp(-x).*sin(x)');[xmax,ymax]=fminbnd(fun1,0,8);ymax=-ymax;fprintf('o:[xmin=%f,ymin=%f]\n',[xmin,ymin]),fprintf('*:[xmax=%f,ymax=%f]\n',[xmax,ymax]),plot(xmin,ymin,'o',xmax,ymax,'p');

例2對(duì)邊長為3米的正方形鐵板，在四個(gè)角剪去相等的正方形以制成方形無蓋水槽，問如何剪法使水槽的容積最大？先編寫M文件fun2.m如下:functionf=fun2(x)f=-(3-2*x).^2*x;主程序?yàn)閙ain2.m:[x,fval]=fminbnd('fun2',0,1.5);xmax=xfmax=-fval運(yùn)算結(jié)果為:xmax=0.5000,fmax=2.0000.即剪掉的正方形的邊長為0.5米時(shí)水槽的容積最大,最大容積為2立方米.解)2)2

例3對(duì)邊長為3米的正方形鐵板，在四個(gè)角剪去相等的正方形以制成方形無蓋水槽，問如何剪法使水槽的容積最大？先編寫M文件fun2.m如下:functionf=fun2(x)f=-(3-2*x).^2*x;主程序?yàn)閙ain2.m:[x,fval]=fminbnd('fun2',0,1.5);xmax=xfmax=-fval運(yùn)算結(jié)果為:xmax=0.5000,fmax=2.0000.即剪掉的正方形的邊長為0.5米時(shí)水槽的容積最大,最大容積為2立方米.解)2)2命令格式為:（1）x=fminunc（fun,X0

）；或x=fminsearch（fun,X0

）（2）x=fminunc（fun,X0

，options）；或x=fminsearch（fun,X0

，options）（3）[x，fval]=fminunc（...）；或[x，fval]=fminsearch（...）（4）[x，fval，exitflag]=fminunc（...）；或[x，fval，exitflag]=fminsearch（5）[x，fval，exitflag，output]=fminunc（...）；或[x，fval，exitflag，output]=fminsearch（...）3.1.2多元函數(shù)無約束優(yōu)化問題標(biāo)準(zhǔn)型為：minF(X)[3]fminunc為中型優(yōu)化算法的步長一維搜索提供了兩種算法，由options中參數(shù)LineSearchType控制：LineSearchType='quadcubic'(缺省值)，混合的二次和三次多項(xiàng)式插值；LineSearchType='cubicpoly'，三次多項(xiàng)式插使用fminunc和fminsearch可能會(huì)得到局部最優(yōu)解.說明:fminsearch是用單純形法尋優(yōu)；fminunc的算法見以下幾點(diǎn)說明：[1]fminunc為無約束優(yōu)化提供了大型優(yōu)化和中型優(yōu)化算法。由options中的參數(shù)LargeScale控制：LargeScale='on'(默認(rèn)值),使用大型算法LargeScale='off'(默認(rèn)值),使用中型算法[2]fminunc為中型優(yōu)化算法的搜索方向提供了4種算法，由options中的參數(shù)HessUpdate控制：HessUpdate='bfgs'（默認(rèn)值），擬牛頓法的BFGS公式；HessUpdate='dfp'，擬牛頓法的DFP公式；HessUpdate='steepdesc'，最速下降法例3minf(x)=(4x12+2x22+4x1x2+2x2+1)*exp(x1)(1)編寫M-文件fun3.m:functionf=fun3(x)f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);(2)輸入M文件main3.m如下:x0=[-1,1];x=fminunc('fun3',x0);y=fun3(x)(3)運(yùn)行結(jié)果:x=0.5000-1.0000y=1.3029e-10例4產(chǎn)銷量的最佳安排

某廠生產(chǎn)一種產(chǎn)品有甲、乙兩個(gè)牌號(hào)，討論在產(chǎn)銷平衡的情況下如何確定各自的產(chǎn)量，使總利潤最大.所謂產(chǎn)銷平衡指工廠的產(chǎn)量等于市場(chǎng)上的銷量.基本假設(shè)(1)價(jià)格與銷量成線性關(guān)系(2)成本與產(chǎn)量成負(fù)指數(shù)關(guān)系

模型建立

若根據(jù)大量的統(tǒng)計(jì)數(shù)據(jù),求出系數(shù)b1=100,a11=1,a12=0.1,b2=280,a21=0.2,a22=2,r1=30,λ1=0.015,c1=20,r2=100,λ2=0.02,c2=30,則問題轉(zhuǎn)化為無約束優(yōu)化問題：求甲,乙兩個(gè)牌號(hào)的產(chǎn)量x1，x2，使總利潤z最大.為簡化模型,先忽略成本,并令a12=0,a21=0,問題轉(zhuǎn)化為求:z1=(b1-a11x1)x1+(b2-a22x2)x2

的極值.顯然其解為x1=b1/(2a11)

=50,x2=b2/(2a22)

=70,令它作為原問題的初始值。總利潤為：z(x1,x2)=(p1-q1)x1+(p2-q2)x2

模型求解(1)建立M-文件fun.m:functionf=fun(x)y1=((100-x(1)-0.1*x(2))-(30*exp(-0.015*x(1))+20))*x(1);y2=((280-0.2*x(1)-2*x(2))-(100*exp(-0.02*x(2))+30))*x(2);f=-y1-y2;(2)輸入命令:x0=[50,70];x=fminunc('fun',x0),z=fun(x)(3)計(jì)算結(jié)果:x=23.9025,62.4977,z=6.4135e+003

即甲的產(chǎn)量為23.9025,乙的產(chǎn)量為62.4977,最大利潤為6413.5。3.2非線性規(guī)劃

3.2.2一般非線性規(guī)劃3.2.1Matlab計(jì)算:二次規(guī)劃用MATLAB軟件求解,其輸入格式如下:（1）x=quadprog(H,c,A,b);（2）x=quadprog(H,c,A,b,Aeq,beq);（3）x=quadprog(H,c,A,b,Aeq,beq,VLB,VUB);（4）x=quadprog(H,c,A,b,Aeq,beq,VLB,VUB,X0);（5）x=quadprog(H,c,A,b,Aeq,beq,VLB,VUB,X0,options);（6）[x,fval]=quaprog(...);（7）[x,fval,exitflag]=quaprog(...);（8）[x,fval,exitflag,output]=quaprog(...);3.2.1Matlab計(jì)算:二次規(guī)劃例1minf(x1,x2)=-2x1-6x2+x12-2x1x2+2x22s.t.x1+x2≤2-x1+2x2≤2x1≥0,x2≥01）寫成標(biāo)準(zhǔn)形式：2）輸入命令：

H=[1-1;-12];c=[-2;-6];A=[11;-12];b=[2;2];Aeq=[];beq=[];VLB=[0;0];VUB=[];[x,z]=quadprog(H,c,A,b,Aeq,beq,VLB,VUB)3）運(yùn)算結(jié)果為：

x=0.66671.3333z=-8.22221）首先建立目標(biāo)函數(shù)f(X)的m文件objFun.m,定義:functionfun=objFun(X);fun=f(X);3.2.2一般非線性規(guī)劃其中X為n維變?cè)蛄浚珿(X)與Ceq(X)均為非線性函數(shù)組成的向量，其它變量的含義與線性規(guī)劃、二次規(guī)劃中相同.用Matlab求解上述問題，基本步驟分三步：3）建立主程序.非線性規(guī)劃求解的函數(shù)是fmincon,命令的基本格式如下：

(1)x=fmincon('fun',X0,A,b)

(2)x=fmincon('fun',X0,A,b,Aeq,beq)

(3)x=fmincon('fun',X0,A,b,Aeq,beq,VLB,VUB)

(4)x=fmincon('fun',X0,A,b,Aeq,beq,VLB,VUB,'nonlcon')(5)x=fmincon('fun',X0,A,b,Aeq,beq,VLB,VUB,'nonlcon',options)

(6)[x,fval]=fmincon(...)

(7)[x,fval,exitflag]=fmincon(...)

(8)[x,fval,exitflag,output]=fmincon(...)輸出極值點(diǎn)M文件迭代的初值參數(shù)說明變量上下限注意：[1]fmincon函數(shù)提供了大型優(yōu)化算法和中型優(yōu)化算法。默認(rèn)時(shí)，若在fun函數(shù)中提供了梯度（options參數(shù)的GradObj設(shè)置為'on'），并且只有上下界存在或只有等式約束，fmincon函數(shù)將選擇大型算法。當(dāng)既有等式約束又有梯度約束時(shí)，使用中型算法。[2]fmincon函數(shù)的中型算法使用的是序列二次規(guī)劃法。在每一步迭代中求解二次規(guī)劃子問題，并用BFGS法更新拉格朗日Hessian矩陣。[3]fmincon函數(shù)可能會(huì)給出局部最優(yōu)解，這與初值X0的選取有關(guān)。1）寫成標(biāo)準(zhǔn)形式：

s.t.

2x1+3x26s.tx1+4x25x1,x20例22）先建立M-文件objFun2.m:functionf=objFun2(x);f=-x(1)-2*x(2)+(1/2)*x(1)^2+(1/2)*x(2)^23）再建立主程序main2.m：

x0=[1;1];A=[23;14];b=[6;5];Aeq=[];beq=[];VLB=[0;0];VUB=[];[x,fval]=fmincon('objFun2',x0,A,b,Aeq,beq,VLB,VUB)4）運(yùn)算結(jié)果為：

x=0.76471.0588fval=-2.02941）先建立M文件objFun3.m,定義目標(biāo)函數(shù):functionf=objFun3(x);f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);

x1+x2=0s.t.1.5+x1x2-x1-x2≤

0-x1x2–10

≤

0例32）再建立M文件mycon.m定義非線性約束：

function[g,ceq]=mycon(x)g=[1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];ceq=[];3）主程序main3.m為:x0=[-1;1];A=[];b=[];Aeq=[11];beq=[0];vlb=[];vub=[];[x,fval]=fmincon('objFun3',x0,A,b,Aeq,beq,vlb,vub,'mycon')4）運(yùn)算結(jié)果為：

x=-1.22501.2250fval=1.8951例4

1）先建立M-文件funGoal.m定義目標(biāo)函數(shù):functionf=funGoal(x);f=-2*x(1)-x(2);2）再建立M文件mycon4.m定義非線性約束：

function[g,ceq]=mycon4(x)g=[x(1)^2+x(2)^2-25;x(1)^2-x(2)^2-7];

3）主程序maincon.m為:x0=[3;2.5];VLB=[00];VUB=[510];[x,fval,exitflag,output]…=fmincon('funGoal',x0,[],[],[],[],VLB,VUB,'mycon4')4）運(yùn)算結(jié)果為:

x=4.00003.0000fval=-11.0000exitflag=1output=iterations:4funcCount:17stepsize:1algorithm:[1x44char]firstorderopt:[]cgiterations:[]例5、應(yīng)用實(shí)例：商品最優(yōu)存儲(chǔ)問題

有一公司，為節(jié)約成本、減少開支，希望盡量減少商品庫存空間，同時(shí)強(qiáng)調(diào)庫存的商品能夠滿足客戶的需求。假設(shè)公司銷售甲乙兩種商品，各種信息符號(hào)如下表所示。

各種符號(hào)所表示的意義xi第i種商品的存儲(chǔ)量ai第i種商品的價(jià)格，a1=9，a2=4bi第i種商品的供給率，b1=3，b2=5hi第i種商品的每單位的存儲(chǔ)費(fèi)用，h1=0.5，h2=0.2ti第i種商品的每單位的存儲(chǔ)空間，t1=2，t2=4T最大存儲(chǔ)空間T=24例5、應(yīng)用實(shí)例：商品最優(yōu)存儲(chǔ)問題根據(jù)歷史數(shù)據(jù)得到商品i的成本可以表示為：問題的目標(biāo)函數(shù)為兩種商品的總成本：第i種商品的存儲(chǔ)空間為：兩種商品所占有的總的存儲(chǔ)空間最大值為T，表達(dá)成約束的形式為：建立非線性規(guī)劃模型為：例5、應(yīng)用實(shí)例：商品最優(yōu)存儲(chǔ)問題代入?yún)?shù)后，非線性規(guī)劃模型為：定義最優(yōu)化問題的目標(biāo)函數(shù)con1_objFun.m如下：functionf=con1_objFun(x)f=27/x(1)+0.25*x(1)+20/x(2)+0.10*x(2);定義約束函數(shù)con1_confun.m如下：function[C,Ceq]=con1_confun(x)C=[2*x(1)+4*x(2)-24];Ceq=[];求解該最優(yōu)化問題的主程序腳本文件main_con1：x0=[3;3];vlb=[0;0];vub=[];options=optimset('algorithm','interior-point');[x,fval,flag,outp]=fmincon('con1_objFun',x0,[],[],[],[],vlb,vub,'con1_confun',options)例5、應(yīng)用實(shí)例：商品最優(yōu)存儲(chǔ)問題運(yùn)行結(jié)果為：x=[5.0968,3.4516]fval=12.7112exitflag=1

output=iterations:7funcCount:24constrviolation:0stepsize:5.8305e-06algorithm:'interior-point'firstorderopt:5.6765e-07cgiterations:0message:[1x777char]例6、應(yīng)用實(shí)例：供應(yīng)與選址某公司有6個(gè)建筑工地要開工，每個(gè)工地的位置（用平面坐標(biāo)系a，b表示，距離單位：千米）及水泥日用量d(噸)由下表給出。目前有兩個(gè)臨時(shí)料場(chǎng)位于A(5,1)，B(2,7)，日儲(chǔ)量各有20噸。假設(shè)從料場(chǎng)到工地之間均有直線道路相連。（1）試制定每天的供應(yīng)計(jì)劃，即從A，B兩料場(chǎng)分別向各工地運(yùn)送多少噸水泥，使總的噸千米數(shù)最小。（2）為了進(jìn)一步減少噸千米數(shù)，打算舍棄兩個(gè)臨時(shí)料場(chǎng)，改建兩個(gè)新的，日儲(chǔ)量各為20噸，問應(yīng)建在何處，節(jié)省的噸千米數(shù)有多大？（1）建立模型

記工地的位置為(ai，bi)，水泥日用量為di，i=1,…,6;料場(chǎng)位置為(xj，yj)，日儲(chǔ)量為ej，j=1,2；從料場(chǎng)j向工地i的運(yùn)送量為Xij。要同時(shí)確定料場(chǎng)的位置(xj,yj)和運(yùn)送量Xij，在同樣條件下使總噸千米數(shù)最小。這是非線性規(guī)劃問題。非線性規(guī)劃模型為：決策變量為：Xij，xj，yj，i=1~6，j=1~2。①先編寫M文件定義目標(biāo)函數(shù)。functionf=con2_objFun(x)a=[1.258.750.55.7537.251.250.754.7556.57.25];

f1=0;fori=1:6s(i)=sqrt((x(13)-a(1,i))^2+(x(14)-a(2,i))^2);

f1=s(i)*x(i)+f1;end

fori=7:12

s(i)=sqrt((x(15)-a(1,i-6))^2+(x(16)-a(2,i-6))^2);f1=s(i)*x(i)+f1;end

f=f1;（2）求解編程X11=x1,X21=x2,X31=x3,X41=x4,X51=x5,X61=x6,X12=x7,X22=x8,X32=x9,X42=x10,X52=x11,,X62=x12x1=X13,y1=X14,x2=X15,y2=X16

②

編寫主程序x0=[35070100406105127]';

A=[11111100000000000000001111110000];e=[20;20];Aeq=[100000100000000001000001000000000010000010000000000100000100000000001000001000000000010000010000];d=[3547611];vlb=[zeros(12,1);-inf;-inf;-inf;-inf];vub=[];[x,fval,exitflag]=fmincon('con2_objFun',x0,A,e,Aeq,d,vlb,vub)(3)計(jì)算結(jié)果為：x=[35471000005115.69464.92677.257.25]fval=89.3118exitflag=1工地123456坐標(biāo)料場(chǎng)A354710(5.6946,4.9267)料場(chǎng)B0000511(7.25,7.25)兩個(gè)料場(chǎng)的位置坐標(biāo)及其供應(yīng)6個(gè)工地的物料量優(yōu)化結(jié)果見下表，最小工作量為89.3。兩個(gè)料場(chǎng)供應(yīng)工地的用量及位置坐標(biāo)4、動(dòng)態(tài)規(guī)劃4.1多階段決策過程最優(yōu)化問題舉例4.2動(dòng)態(tài)規(guī)劃的基本概念與建模編程4.3動(dòng)態(tài)規(guī)劃案例分析動(dòng)態(tài)規(guī)劃是解決多階段決策過程的最優(yōu)化問題的一種方法。

多階段決策過程的決策問題:可以按時(shí)間、空間等標(biāo)識(shí)把它分為很多階段，每一階段都需作出決策，使得整個(gè)過程達(dá)到最優(yōu)。

動(dòng)態(tài)規(guī)劃方法把多階段決策問題變成一系列互相聯(lián)系的單階段問題，解決了一系列較容易的單階段問題，也就解決了困難的多階段決策問題。例1

商品定價(jià)問題

某廠要確定一種新產(chǎn)品在今后五年內(nèi)的價(jià)格，并擬訂只在5，6，7，8元這四種單價(jià)中選擇。根據(jù)預(yù)測(cè)，今后五年內(nèi)不同價(jià)格下產(chǎn)品的每年盈利額(萬元)如下表。如果需要調(diào)整價(jià)格，各相鄰年度的價(jià)格調(diào)整(增減)不得超過1元。問今后五年內(nèi)產(chǎn)品每年定價(jià)各為多少，可使五年的總利潤最大？4.1多階段決策過程最優(yōu)化問題舉例利潤年價(jià)格

123455

6

7

892458

75864

65973

87664解：用標(biāo)號(hào)法求解。將問題分為五個(gè)階段，逆序遞推，求利潤最大，計(jì)算結(jié)果如下圖所示。年價(jià)格

123455

6

7

892458

75864

65973

87664[8][13][18][24][37][4][14][22][28][35][4][10][17][30][38][3][11][23][28][36]4.2動(dòng)態(tài)規(guī)劃的基本概念與建模編程4.2.1基本概念（1）階段k根據(jù)決策問題性質(zhì)，可將決策過程（按空間位置、時(shí)間進(jìn)程、工序）恰當(dāng)劃分為若干相聯(lián)系的階段，總決策是各階段的序列決策之合。各階段用k表示，稱為階段變量。如：例1劃分成4個(gè)階段，CiDj位于第三階段例2劃分成5個(gè)階段。（2）決策xk

從某一狀態(tài)向下一狀態(tài)過度時(shí)所做的選擇，記為xk，xk

稱為決策變量。

例1：x2(B1)=C1，另外，還可以有x2(B1)=C2，x2(B1)=C3。在狀態(tài)sk下，允許采取決策的全體稱為決策允許集合，記為Dk(sk)

。如

D2(s2)=D2(B3)={B3C1,B3C2,B3C3}。（3）狀態(tài)sk階段開始的狀況或條件稱為狀態(tài)，是表示決策過程當(dāng)前特征的量，是決策的出發(fā)點(diǎn)。表示狀態(tài)的變量(sk)，稱為狀態(tài)變量，其應(yīng)具有無后效性。狀態(tài)的描述可以是數(shù)量，也可以是字符，數(shù)量狀態(tài)可以是連續(xù)的，也可以是離散的。

例1中：s2=B1表示第2階段當(dāng)前在B1

位置，s2還可以是另外兩個(gè)狀態(tài)B2和B3。例2中：每年的初始價(jià)格sk

為狀態(tài)變量。（4）狀態(tài)轉(zhuǎn)移方程sk+1=Tk(sk，xk)下一階段的狀態(tài)與本階段狀態(tài)及決策的函數(shù)關(guān)系。多階段決策過程中，給定第k階段的狀態(tài)sk，以及本階段的決策xk后，第k+1階段的狀態(tài)sk+1也隨之確定：sk+1=Tk(sk,xk)例如，s3=C1=T2(B2，C1)狀態(tài)

s1階段1T1決策x1狀態(tài)

s2決策x2階段2T2狀態(tài)

s3...狀態(tài)

sk決策xk階段kTk狀態(tài)

sk+1...狀態(tài)

sn決策xn階段nTn狀態(tài)

sn+1（5）多階段決策過程：k-n（后部）子過程全過程（6）策略Pk,n(sk)從第k階段開始到最后第n階段的決策序列，稱k子策略。P1,n(s1)

稱為全過程策略。例如，P1,5(A)=

AB2C1D1E，就是一個(gè)全過程策略，而且是最優(yōu)策略；而P3,5(C2)=C2D1E是一個(gè)3-5

子策略。

（7）階段指標(biāo)函數(shù)rk(sk，xk)從狀態(tài)sk出發(fā)，由決策xk所產(chǎn)生的第k階段效益稱為第k階段指標(biāo)，這種關(guān)系稱為階段指標(biāo)函數(shù)。例如r2(B3，C2)=12。（8）過程指標(biāo)函數(shù)Rk,n(sk，xk，xk+1，…，xn)從狀態(tài)sk出發(fā)，選擇決策xk，xk+1，…，xn所產(chǎn)生的過程（效益）指標(biāo)的關(guān)系。例如R2,5(B3，C2，D1，E

)=23。動(dòng)態(tài)規(guī)劃要求過程指標(biāo)函數(shù)具有可分離性。即Rk，n(sk，xk，xk+1，…，xn)=rk(sk，xk)+

Rk+1(sk+1，xk+1，…，xn)，稱指標(biāo)具有可加性，或Rk，n(sk，xk，xk+1，…，xn)=rk(sk，xk)×Rk+1(sk+1，xk+1，…，xn)，稱指標(biāo)具有可乘性。（9）最優(yōu)指標(biāo)函數(shù)fk(sk)

最優(yōu)指標(biāo)函數(shù)fk(sk)：從狀態(tài)sk出發(fā)，對(duì)所有的策略Pk,n，過程指標(biāo)Rk,n的最優(yōu)值，即對(duì)于可加性指標(biāo)函數(shù)，最優(yōu)指標(biāo)函數(shù)fk(sk)為對(duì)于可乘性指標(biāo)函數(shù)，最優(yōu)指標(biāo)函數(shù)fk(sk)為上兩式中“opt”表示“max”或“min”。上式稱為動(dòng)態(tài)規(guī)劃最優(yōu)指標(biāo)的遞推方程，是動(dòng)態(tài)規(guī)劃的基本方程。終端條件：為了使以上的遞推方程有遞推的起點(diǎn)，必須要設(shè)定最優(yōu)指標(biāo)的終端條件，一般對(duì)于可加性指標(biāo)設(shè)n+1階段的最優(yōu)指標(biāo)fn+1(sn+1)=0；對(duì)于可乘性指標(biāo)設(shè)n+1階段的最優(yōu)指標(biāo)fn+1(sn+1)=1。4.2.2動(dòng)態(tài)規(guī)劃的數(shù)學(xué)模型構(gòu)建步驟1)階段變量：k=1,…,n,n+1(n+1為邊界條件設(shè)定)2)第k階段的決策變量為xk，允許決策集合為Dk(sk)3)第k階段的狀態(tài)變量為sk，狀態(tài)集合為Sk狀態(tài)變量應(yīng)具備的特點(diǎn)：①能夠描述決策過程的演變特性，反映決策的依據(jù)及與決策結(jié)果的聯(lián)系；②滿足無后效性，不跨階段影響決策。即某階段狀態(tài)確定后，未來的發(fā)展不受以前各階段狀態(tài)的影響；③可知性，即能夠獲得狀態(tài)變量的值。4)狀態(tài)轉(zhuǎn)移方程：sk+1=T（sk，xk）5)階段效益(指標(biāo)函數(shù))rk=rk(sk，xk)6)過程指標(biāo)函數(shù)Rk，n(sk，xk，xk+1，…，xn)=∑ri(si，xi)7)第k階段的最優(yōu)指標(biāo)函數(shù)fk=fk(sk)則基本方程為：

fk(sk)=Opt{Rk,n(sk,Pk,n(xk))}Pk,n=Opt{rk(sk,xk)+Rk+1,n(sk+1,Pk+1,n)}(xk,Pk+1,n)=Opt{rk(sk,xk)+OptRk+1,n(sk+1,Pk+1,n)}xkPk+1,n=Opt{rk(sk,xk)+fk+1,n(sk+1)}k=n,n-1,…,2,1xkfn+1=0(邊界條件)ni=k4.2.3動(dòng)態(tài)規(guī)劃的MATLAB程序function[p_opt,fval]=zxf_dynprog(s,DecisFun,SubObjfun,TransFun,Objfun)%輸出參數(shù)：%p_opt——由4列構(gòu)成的最優(yōu)策略信息，p_opt=[序號(hào),狀態(tài),最優(yōu)策略,指標(biāo)函數(shù)值]%fval——向量各元素分別表示p_opt各最優(yōu)策略組對(duì)應(yīng)始端狀態(tài)s的最優(yōu)函數(shù)值%輸入?yún)?shù)：%s——狀態(tài)(m)*階段(n)矩陣，一列代表一個(gè)階段的所有狀態(tài)%DecisFun(k,s)——定義由階段k、狀態(tài)s確定的允許決策（集合）向量%TransFun(k,s,x)——狀態(tài)轉(zhuǎn)移函數(shù)sk+1=T(sk,xk)，s是階段狀態(tài)，x是決策變量%SubObjfun(k,s,x)——定義階段k、狀態(tài)s、決策x對(duì)應(yīng)的階段效益指標(biāo)函數(shù)r%Objfun(v,f)——第k-n過程指標(biāo)函數(shù)Rk,n，當(dāng)Rk,n=Objfun(v,f)=v+f時(shí)可省略

[m,n]=size(s);%計(jì)算狀態(tài)*階段的矩陣尺寸s_is_s=~isnan(s);%狀態(tài)*階段矩陣的非空矩陣f_opt=nan*ones(m,n);%對(duì)應(yīng)狀態(tài)*階段矩陣的最優(yōu)過程指標(biāo)矩陣d_opt=f_opt;%對(duì)應(yīng)狀態(tài)*階段矩陣的決策矩陣

%先計(jì)算最后階段的相關(guān)值Index_sn=find(s_is_s(:,n));%取出最后階段狀態(tài)向量索引num_sn=length(Index_sn);%計(jì)算最后階段狀態(tài)個(gè)數(shù)fori=1:num_sn%執(zhí)行DecisFun，獲取n階段s(i,n)狀態(tài)的決策(集合)向量x=feval(DecisFun,n,s(Index_sn(i),n));num_x=length(x);%計(jì)算決策向量元素個(gè)數(shù)R_value=-inf;forj=1:num_x%循環(huán)求出當(dāng)前狀態(tài)下效益指標(biāo)最大的決策%計(jì)算在Index_sn(i)的狀態(tài)下x(j)決策的指標(biāo)值f_tmp=feval(SubObjfun,n,s(Index_sn(i),n),x(j));iff_tmp>R_value%當(dāng)指標(biāo)值大時(shí)f_opt(Index_sn(i),n)=f_tmp;%保存指標(biāo)值d_opt(Index_sn(i),n)=x(j);%保存決策值R_value=f_tmp;%保存指標(biāo)值end;end;end

%逆推計(jì)算(n-1)~1各階段指標(biāo)的遞歸調(diào)用程序fork=n-1:-1:1Index_sk=find(s_is_s(:,k));%取出第k階段狀態(tài)向量索引num_sk=length(Index_sk);%計(jì)算第k階段狀態(tài)個(gè)數(shù)fori=1:num_sk%執(zhí)行DecisFun，獲取k階段s(i,k)狀態(tài)的決策（集合）向量x=feval(DecisFun,k,s(Index_sk(i),k));num_xk=length(x);%計(jì)算決策向量元素個(gè)數(shù)R_value=-inf;forj=1:num_xk%循環(huán)求出當(dāng)前狀態(tài)下效益指標(biāo)最大的決策%計(jì)算在Index_sk(i)狀態(tài)下x(j)決策的階段指標(biāo)r_value=feval(SubObjfun,k,s(Index_sk(i),k),x(j));%計(jì)算在Index_sk(i)狀態(tài)下x(j)決策的狀態(tài)轉(zhuǎn)移s_next=feval(TransFun,k,s(Index_sk(i),k),x(j));tmp0=s(:,k+1)-s_next;index_sk1=find(tmp0==0);%獲取狀態(tài)對(duì)應(yīng)的索引if~isempty(index_sk1),ifnargin<5f_tmp=r_value+f_opt(index_sk1(1),k+1);%計(jì)算過程指標(biāo)（默認(rèn)表達(dá)式）elsef_tmp=feval(Objfun,r_value,f_opt(index_sk1(1),k+1));%計(jì)算過程指標(biāo)endiff_tmp>R_value%當(dāng)過程指標(biāo)值大時(shí)f_opt(Index_sk(i),k)=f_tmp;%保存過程指標(biāo)值d_opt(Index_sk(i),k)=x(j);%保存決策值R_value=f_tmp;%保存過程指標(biāo)值end;end;end;end;end;fval=f_opt(Index_sn,1);fval=fval(find(~isnan(fval)),1);%輸出fval4.2.3動(dòng)態(tài)規(guī)劃的MATLAB程序%記錄最優(yōu)決策、最優(yōu)軌線和相應(yīng)指標(biāo)函數(shù)值p_opt=[];tmp_s=[];tmp_d=[];tmp_r=[];Index_s1=find(s_is_s(:,1));%第1階段有效狀態(tài)索引num_s1=length(Index_s1);%第1階段有效狀態(tài)數(shù)目fori=1:num_s1,tmp_d(i)=d_opt(Index_s1(i),1);%取第1階段i狀態(tài)下的決策tmp_s(i)=s(Index_s1(i),1);%取第1階段i狀