本文格式為Word版，下載可任意編輯——單純形法線性規(guī)劃實(shí)踐報告..

...

線性規(guī)劃——單純形法程序設(shè)計

1.試驗?zāi)康模?/p>

(1)使學(xué)生在程序設(shè)計方面得到進(jìn)一步的訓(xùn)練;，把握Matlab(C或VB)語言進(jìn)行程序設(shè)計中一些常用方法。

(2)使學(xué)生對線性規(guī)劃的單純形法有更深的理解.

2.問題述

本試驗主要編寫一般線性規(guī)劃問題的計算程序:

Min

s.t.

x

引入松弛變量將其化為一般標(biāo)準(zhǔn)型線性規(guī)劃問題：

Min

s.t.Ax=b;

x

A為m*n的矩陣，有m個約束，n個變量。

求解上述線性規(guī)劃采用單純形算法，初始可行基由引入的m個人工變量對應(yīng)的單位陣組成，并采用大M算法

3.算法描述

（1）將引入的人工變量對應(yīng)的單位陣作為初始可行基，則原線性規(guī)劃問題構(gòu)造出下面的新線性規(guī)劃問題：

（2）通過判別數(shù)計算公式可求出n+m個變量的判別數(shù)，若全部判別數(shù),則得到一個最優(yōu)基本可行解，運(yùn)算終止；否則，轉(zhuǎn)到下一步

（3）找出判別數(shù)為負(fù)的最小判別數(shù)，其對應(yīng)的變量為入基變量，記錄下來標(biāo)為k，然后看其對應(yīng)的列向量，若中的所有元，則原線性規(guī)劃無最優(yōu)解，否則，轉(zhuǎn)下一步

（4）計算，則為離基變量，然后對A進(jìn)行初等變換，計算

（5）用入基變量與出基變量對應(yīng)的列向量、判別、對應(yīng)的系數(shù)均對換，則可用計算機(jī)編程循環(huán)以上步驟直至得出結(jié)果

4.計算程序（matlab）

程序保存為linpro.m文件

5.算例驗證

在窗口輸入：

Aeq=[1-1-1-110;0-1-1101;111100];

b=[0;0;1];

c0=[-0.15-0.1-0.08-0.1200];

linpro(Aeq,b,c0)

10000

10000

-0.1300

說明通過三次迭代找到最優(yōu)解為-0.13.

用Matlab求解線性規(guī)劃的命令linprog的計算結(jié)果：

f=[-0.15;-0.1;-0.08;-0.12];

A=[1-1-1-1

0-1-11];

b=[0;0];

Aeq=[1111];

beq=[1];

lb=zeros(4,1);

然后調(diào)用linprog函數(shù)：

[x,fval]=linprog(f,A,b,Aeq,beq,lb)；

x=

0.5000

0.2500

0.0000

0.2500

fval=

-0.1300

最優(yōu)值為-0.13，與上面的結(jié)果一致，說明程序正確。

單變量在單峰區(qū)間上的微小點(diǎn)——黃金分割法

問題述

設(shè)f(x)為區(qū)間[a,b]上的下單峰函數(shù)，在區(qū)間取兩點(diǎn)x1,x2，使每次迭代都把區(qū)間縮短率定為0.618.x1x2，若f(x1)f(x2),則[a,x2];若f(x1)f(x2),則.這樣區(qū)間就減小了，直至找到最優(yōu)解。

算法描述

令x2=a+0.618(b-a),f2=f(x2);

令x1=a+0.382(b-a),f1=f(x1);

若b-a,則找到最優(yōu)解=（a+b）/2；否則轉(zhuǎn)下一步；

若f1f2,則b=x2,x2=x1,f2=f1,轉(zhuǎn)其次步；

若f1=f2,則a=x1,b=x2,轉(zhuǎn)第一步；

若f1f2,則a=x1,x1=x2,f1=f2,轉(zhuǎn)第三步；

（5）令x2=a+0.618(b-a),f2=f(x2),轉(zhuǎn)第三步。

3.程序代碼（C++）

以函數(shù)f(x)=為例編寫：

#includeiostream

usingnamespacestd;

doublea,b,t;

doublef(doublex)

{

doubled;

d=x*x-x+2;

returnd;}

doubleminf(doublea,doubleb)

{

doublex1,x2,f1,f2;

x1=a+0.382*(b-a);

x2=a+0.618*(b-a);

f1=f(x1);

f2=f(x2);

while(b-a1e-4)

{

if(f1f2)

{

b=x2;x2=x1;f2=f1;

x1=a+0.382*(b-a);

f1=f(x1);}

elseif(f1f2)

{

a=x1;x1=x2;f1=f2;

x2=a+0.618*(b-a);

f2=f(x2);}

else

{a=x1;b=x2;

x1=a+0.382*(b-a);

x2=a+0.618*(b-a);

f1=f(x1);

f2=f(x2);}

t=(a+b)/2.0;

}

returnt;

}

voidmain()

{

doublek;

cout請輸入單峰區(qū)間的端點(diǎn)a和b：;

cinab;

k=minf(a,b);

cout最優(yōu)值為：k最優(yōu)解為：f(k)endl;

}

算例驗證

函數(shù)f(x)=在區(qū)間[-1,3]上的最小值，程序運(yùn)算如下：

在matlab命令窗口輸入如下：

fun=x^2-x+2;

[x,fval]=fminbnd(fun,-1,3)

運(yùn)算結(jié)果如下：

x=0.5000

fval=1.7500

兩者運(yùn)行結(jié)果完全一致，說明程序正確。

三、運(yùn)用非線性規(guī)劃建模的實(shí)例

1.問題描述：

試設(shè)計一壓縮圓柱螺旋彈簧，要求其質(zhì)量最小。彈簧材料為65Mn，最大工作載荷Pmax=40N，最小工作載荷為0，載荷變化頻率fr=25Hz，彈簧壽命為104h，彈簧鋼絲直徑d的取值圍為1-4mm，中徑D2的取值圍為10-30mm，工作圈數(shù)n不應(yīng)小于4.5圈，彈簧旋繞比C不應(yīng)小于4，彈簧一端固定，一端自由，工作溫度為50C，彈簧變形量不小于10mm。

2.模型建立

此題的優(yōu)化目標(biāo)是使彈簧質(zhì)量最小，圓柱螺旋彈簧的質(zhì)量可以表示為

式中，--彈簧材料的密度，對于鋼材=7.8×10-6kg/mm3；

n—工作圈數(shù)；

n2—死圈數(shù)，常取n2=1.5-2.5，現(xiàn)取n2=2；

D2—彈簧中徑，mm；

d—彈簧鋼絲直徑，mm。

將已知參數(shù)代入公式，進(jìn)行整理以后得到問題的目標(biāo)函數(shù)為

根據(jù)彈簧性能和結(jié)構(gòu)上的要求，可寫出問題的約束條件：

強(qiáng)度條件

剛度條件

穩(wěn)定性條件

不發(fā)生共振現(xiàn)象，要求

彈簧旋繞比的限制

對d,n,D2的限制

且應(yīng)取標(biāo)準(zhǔn)值，即1.0,1.2,1.6,2.0,2.5,3.0,3.5,4.0mm等。

由上可知，該壓縮圓柱螺旋彈簧的優(yōu)化設(shè)計是一個三維的約束優(yōu)化問題，其數(shù)學(xué)模型為：

取初始設(shè)計參數(shù)為X(0)=[2.0,5.0,25.0]T

首先編寫目標(biāo)函數(shù)的M文件opt25_3.m，返回x處的函數(shù)值f。

functionf=myfun(x)

f=0.192457*1e-4*(x(2)+2)*x(1)^2*x(3);

由于約束條件中有非線性約束，所以需要編寫一個描述非線性約束條件的M文件opt25_3c.m：

function[c,ceq]=mycon(x)

c(1)=350-163*x(1)^(-2.86)*x(3)^0.86;

c(2)=10-0.4*0.01*x(1)^(-4)*x(2)*x(3)^3;

c(3)=(x(2)+1.5)*x(1)+0.44*0.01*x(1)^(-4)*x(2)*x(3)^3-3.7*x(3);

c(4)=375-0.356*1e6*x(1)*x(2)^(-1)*x(3)^(-2);

c(5)=4-x(3)/x(1);

然后設(shè)置線性約束的系數(shù)：

A=[-100

100

0–10

010

00–1

001];

b=[-1;4;-4.5;50;-10;30];

下一步給定初值，給定變量的下限約束，并調(diào)用優(yōu)化過程(磁盤中M文件為opt25_3.m)

x0=[2.0;5.0;25.0];

lb=zeros(3,1);

[x,fval,exitflag,output,lambda]=fmincon(opt25_3o,x0,A,b,[],[],lb,[],opt25_3c)

計算結(jié)果為：

x=

1.6564

4.5000

16.1141

fval=

0.0055

exitflag=

1

output=

iterations:3

funcCount:16

stepsize:1

algorithm:medium-scale:SQP,Quasi-Newton,li