. .
. . .
线性规划——单纯形法程序设计
1.实验目的：
(1)使学生在程序设计方面得到进一步的训练;，掌握Matlab (C或VB)语言进行程序设计中一些常用方法。
(2)使学生对线性规划的单纯形法有更深的理解.
2.问题述
本实验主要编写一般线性规划问题的计算程序:
Min
s.t.
x
引入松弛变量将其化为一般标准型线性规划问题：
Min
s.t. Ax=b;
x
A为m*n的矩阵，有m个约束，n个变量。
求解上述线性规划采用单纯形算法，初始可行基由引入的m个人工变量对应的单位阵组成，并采用大M算法
3.算法描述
（1）将引入的人工变量对应的单位阵作为初始可行基，则原线性规划问题构造出下面的新线性规划问题：
（2）通过判别数计算公式可求出n+m个变量的判别数，若全部判别数,则得到一个最优基本可行解，运算结束；否则，转到下一步
（3）找出判别数为负的最小判别数，其对应的变量为入基变量，记下标为k，然后看其对应的列向量，若中的所有元，则原线性规划无最优解，否则，转下一步
（4）计算，则为离基变量，然后对A进行初等变换，计算
（5）用入基变量与出基变量对应的列向量、判别、对应的系数均对换，则可用计算机编程循环以上步骤直至得出结果
4.计算程序（matlab）
程序保存为linpro.m文件
5.算例验证
在窗口输入：
Aeq=[1 -1 -1 -1 1 0;0 -1 -1 1 0 1;1 1 1 1 0 0];
b=[0;0;1];
c0=[-0.15 -0.1 -0.08 -0.12 0 0];
linpro(Aeq,b,c0)
10000
10000
-0.1300
说明通过三次迭代找到最优解为-0.13.
用Matlab求解线性规划的命令linprog的计算结果：
f = [-0.15;-0.1;-0.08;-0.12];
A = [1 -1 -1 -1
0 -1 -1 1];
b = [0; 0];
Aeq=[1 1 1 1];
beq=[1];
lb = zeros(4,1);
然后调用linprog函数：
[x,fval] = linprog(f,A,b,Aeq,beq,lb)；
x =
0.5000
0.2500
0.0000
0.2500
fval =
-0.1300
最优值为-0.13，与上面的结果一致，说明程序正确。
单变量在单峰区间上的极小点——黄金分割法
问题述
设f(x)为区间[a,b]上的下单峰函数，在区间取两点x1,x2，使每次迭代都把区间缩短率定为0.618.x1<x2，若f(x1)<f(x2),则[a,x2];若f(x1)>f(x2),则.这样区间就减小了，直至找到最优解。
算法描述
令x2=a+0.618(b-a),f2=f(x2);
令x1=a+0.382(b-a),f1=f(x1);
若b-a<,则找到最优解=（a+b）/2；否则转下一步；
若f1<f2,则b=x2,x2=x1,f2=f1,转第二步；
若f1=f2,则a=x1,b=x2,转第一步；
若f1>f2,则a=x1,x1=x2,f1=f2,转第三步；
（5）令x2=a+0.618(b-a),f2=f(x2),转第三步。
3.程序代码（C++）
以函数f(x)=为例编写：
#include<iostream>
using namespace std;
double a,b,t;
double f(double x)
{
double d;
d=x*x-x+2;
return d;}
double minf(double a,double b)
{
double x1,x2,f1,f2;
x1=a+0.382*(b-a);
x2=a+0.618*(b-a);
f1=f(x1);
f2=f(x2);
while(b-a>1e-4)
{
if(f1<f2)
{
b=x2;x2=x1;f2=f1;
x1=a+0.382*(b-a);
f1=f(x1);}
else if(f1>f2)
{
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;
}
return t;
}
void main()
{
double k;
cout<<"请输入单峰区间的端点a和b：";
cin>>a>>b;
k=minf(a,b);
cout<<"最优值为："<<k<<"最优解为："<<f(k)<<endl;
}
算例验证
函数f(x)=在区间[-1,3]上的最小值，程序运算如下：
在matlab命令窗口输入如下：
fun='x^2-x+2';
[x,fval]=fminbnd(fun,-1,3)
运算结果如下：
x =0.5000
fval = 1.7500
两者运行结果完全一致，说明程序正确。
三、运用非线性规划建模的实例
1.问题描述：
试设计一压缩圆柱螺旋弹簧，要求其质量最小。弹簧材料为65Mn，最大工作载荷Pmax=40N，最小工作载荷为0，载荷变化频率fr=25Hz，弹簧寿命为104h，弹簧钢丝直径d的取值围为1-4mm，中径D2的取值围为10-30mm，工作圈数n不应小于4.5圈，弹簧旋绕比C不应小于4，弹簧一端固定，一端自由，工作温度为50C，弹簧变形量不小于10mm。
2.模型建立
本题的优化目标是使弹簧质量最小，圆柱螺旋弹簧的质量可以表示为
式中，--弹簧材料的密度，对于钢材=7.8×10-6kg/mm3；
n—工作圈数；
n2—死圈数，常取n2=1.5-2.5，现取n2=2；
D2—弹簧中径，mm；
d—弹簧钢丝直径，mm。
将已知参数代入公式，进行整理以后得到问题的目标函数为
根据弹簧性能和结构上的要求，可写出问题的约束条件：
强度条件
刚度条件
稳定性条件
不发生共振现象，要求
弹簧旋绕比的限制
对d,n,D2的限制
且应取标准值，即1.0,1.2,1.6,2.0,2.5,3.0,3.5,4.0mm等。
由上可知，该压缩圆柱螺旋弹簧的优化设计是一个三维的约束优化问题，其数学模型为：
取初始设计参数为X(0)=[2.0,5.0,25.0]T
首先编写目标函数的M文件opt25_3.m，返回x处的函数值f。
function f = 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);
然后设置线性约束的系数：
A=[-1 0 0
1 0 0
0 –1 0
0 1 0
0 0 –1
0 0 1];
b=[-1;4;-4.5;50;-10;30];
下一步给定初值，给定变量的下限约束，并调用优化过程(磁盘中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)
计算结果为：
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, line-search'
firstorderopt: []
cgiterations: []
所以当弹簧钢丝的直径d、工作圈数n及中径D2分别取1.6564、4.5000和16.1141时弹簧质量最小，为5.5克。考虑到实际情况，各参数可分别取1.6、5.0和16.0
x=[1.6 5.0 16.0];
z=optim253(x)
z =
0.0055
故此时弹簧质量仍为5.5克。