最优控制大作业.docx
《最优控制大作业.docx》由会员分享，可在线阅读，更多相关《最优控制大作业.docx（34页珍藏版）》请在冰豆网上搜索。
最优控制大作业
单纯行法
第一题和第二题采用了单纯形法进行解决，单纯形法的理论依据是：
线形规划问题的可行域是n维向量空间Rn中的多面凸集，其最优值如果存在必在该凸集的某顶点处达到。
顶点所对应的可行解称为基本可行解。
单纯形法的基本思想是：
先找出一个基本可行解，对它进行鉴别，看是否是最优解；若不是，则按照一定法则转换到另一改进的基本可行解,再鉴别；若仍不是，则再转换,按此重复进行。
因基本可行解的个数有限，故经有限次转换必能得出问题的最优解。
如果问题无最优解也可用此法判别。
单纯形法的一般解题步骤可归纳如下：
1.把线性规划问题的约束方程组表达成典范型方程组，找出基本可行解作为初始基本可行解。
2.若基本可行解不存在，即约束条件有矛盾，则问题无解。
3.若基本可行解存在，从初始基本可行解作为起点，根据最优性条件和可行性条件，引入非基变量取代某一基变量，找出目标函数值更优的另一基本可行解。
4.按步骤3进行迭代,直到对应检验数满足最优性条件（这时目标函数值不能再改善），即得到问题的最优解。
5.若迭代过程中发现问题的目标函数值无界，则终止迭代。
用单纯形法求解线性规划问题所需的迭代次数主要取决于约束条件的个数。
1某工厂生产A和B两种产品。
已知制造A产品，每公斤要用煤9吨、电力4千万、劳力3个；制造产品B，每公斤要用煤4吨、电力5千瓦-劳力10个。
又知制成产品A每公斤的产值是7万元；B每公斤的产值是12万元。
现该厂只有煤360吨、电力200千瓦、劳力300个。
问在这种条件下，应该生产A、B产品各多少才能使产值为最高。
试写出其数学模型，即约束方程和目标函数，并利用单纯形法求解该线性规划问题。
解：
设生产A、B产品各
吨
使
引入附加变量
，使不等式约束变为等式约束
程序清单如下：
#include
#include
usingnamespacestd;
intvarIn（doubledelta[5]）;//计算进基变量
intvarOut（doublesita[3],doublea[3][5],doubleb[3],doubledelta[5]）;//计算出基变量
voidcal（doublesita[3],doublea[3][5],doubleb[3],doubledelta[5]）;//计算方程组系数和判别数
double*A;//产品A的产量
double*B;//产品B的产量
intmain（）
{
doublemax=0;//最高产值
doubledelta[5]={-7,-12,0,0,0};//判别数
doubleb[3]={360,200,300};//基可行解
doublea[3][5]={9,4,1,0,0,//系数矩阵A
4,5,0,1,0,
3,10,0,0,1};
doublesita[3];//出基变量判别数
intin;//进入基变量位置
intout;//出入基变量位置
intflag=0;//判别数是否全部非负的标志
while（!
flag）
{
in=varIn（delta）;
out=varOut（sita,a,b,delta）;
cal（sita,a,b,delta）;
flag=1;
for（inti=0;i<5;i++）
{
if（delta[i]<0）//是否所有的判别数都大于0
flag=0;
}
}
cout<<"生产A和B产品的产量和最大产值：
"<cout<<"*****************************"<max=7*（*A）+12*（*B）;
cout<<"产值最大时A产品的产量:
x="<<*A<<","<cout<<"产值最大时B产品的产量:
y="<<*B<<","<cout<<"最大产值为：
optimal="<cout<<"*****************************"<return0;
}
intvarIn（doubledelta[5]）
{
intk=0;
doublemindelta=delta[0];
for（inti=1;i<5;i++）
{
if（delta[i]{
mindelta=delta[i];
k=i;
}
}
returnk;
}
intvarOut（doublesita[3],doublea[3][5],doubleb[3],doubledelta[5]）
{
intl=0;
doubleminsita=1000;
for（inti=1;i<3;i++）
{
if（a[i][varIn（delta）]>0）
sita[i]=b[i]/a[i][varIn（delta）];
if（sita[i]{
minsita=sita[i];
l=i;
}
}
returnl;
}
voidcal（doublesita[3],doublea[3][5],doubleb[3],doubledelta[5]）
{
inti,j,k,l;
doublea1[3][5];
doubleb1[3];
doubledelta1[5];
k=varIn（delta）;
l=varOut（sita,a,b,delta）;
for（i=0;i<3;i++）
{
b1[i]=b[i];
for（j=0;j<5;j++）
{
a1[i][j]=a[i][j];
delta1[j]=delta[j];
}
}
for（i=0;i<3;i++）
{
for（j=0;j<5;j++）
{
if（i!
=l）
{
a[i][j]=a1[i][j]-a1[l][j]*a1[i][k]/a1[l][k];
b[i]=b1[i]-b1[l]*a1[i][k]/a1[l][k];
}
else
{
a[i][j]=a1[l][j]/a1[l][k];
b[i]=b1[l]/a1[l][k];
}
delta[j]=delta1[j]-a1[l][j]*delta1[k]/a1[l][k];
}
}
if（k==0）
{
A=&b[l];
}
elseif（k==1）
{
B=&b[l];
}
}
程序运行结果为：
2某车间有一批长度为180公分的钢管（数量充分多）。
为制造零件的需要，要将其截成3种不同长度的管料：
70公分、52公分、35公分。
经过试算，可以由以下8种截法：
（一）
（二）
（三）
（四）
（五）
（六）
（七）
（八）
70公分
2
1
1
1
0
0
0
0
52公分
0
2
1
0
3
2
1
0
35公分
1
0
1
3
0
2
3
5
边料
5cm
6cm
23cm
5cm
24cm
6cm
23cm
5cm
问应如何同时采取若干种截法配合起来，既满足下列情况的配套要求，而又使总的边料为最少。
（a）要求：
70公分长的管料不少于100根，52公分长的管料不少于100根，35公分长的管料不多于100根。
（b）要求：
70公分、52公分、35公分长的管料均不少于100根。
（c）要求：
70公分长的管料不少于100根，52公分长的管料不多于150根，35公分长的管料不少于100根。
试分别写出上述3种要求下的约束条件和目标函数。
利用单纯形法由计算机分别求解3个线性规划问题。
解：
（a）设八种截法各截出的管料数为：
使
引入附加变量
和人为变量
使不等式约束变为等式约束
（b）设八种截法各截出的管料数为：
使
引入附加变量
和人为变量
使不等式约束变为等式约束
使
（c）设八种截法各截出的管料数为：
使
引入附加变量
和人为变量
，使不等式约束变为等式约束
使
程序清单如下：
#include
#include
#include
doubleM=10000000;
staticdoubleA[20][20];
doublefunction（doublex1,doublex2,doublex3,doublex4,doublex5,
doublex6,doublex7,doublex8）
{
doubleval;
val=5*x1+6*x2+23*x3+5*x4+24*x5+6*x6+23*x7+5*x8;
returnval;
}
voidinit_Array（intsign）
{
switch（sign）
{
case1:
A[0][1]=5;A[0][2]=6;A[0][3]=23;A[0][4]=5;A[0][5]=24;
A[0][6]=6;A[0][7]=23;A[0][8]=5;
A[1][0]=100;A[1][1]=2;A[1][2]=1;A[1][3]=1;A[1][4]=1;
A[1][5]=0;A[1][6]=0;A[1][7]=0;A[1][8]=0;A[1][9]=-1;
A[1][10]=0;A[1][11]=0;
A[2][0]=100;A[2][1]=0;A[2][2]=2;A[2][3]=1;A[2][4]=0;
A[2][5]=3;A[2][6]=2;A[2][7]=1;A[2][8]=0;A[2][9]=0;
A[2][10]=-1;A[2][11]=0;
A[3][0]=100;A[3][1]=1;A[3][2]=0;A[3][3]=1;A[3][4]=3;
A[3][5]=0;A[3][6]=2;A[3][7]=3;A[3][8]=5;A[3][9]=0;
A[3][10]=0;A[3][11]=1;
break;
case2:
A[0][1]=5;A[0][2]=6;A[0][3]=23;A[0][4]=5;A[0][5]=24;
A[0][6]=6;A[0][7]=23;A[0][8]=5;
A[1][0]=100;A[1][1]=2;A[1][2]=1;A[1][3]=1;A[1][4]=1;
A[1][5]=0;A[1][6]=0;A[1][7]=0;A[1][8]=0;A[1][9]=-1;
A[1][10]=0;A[1][11]=0;
A[2][0]=100;A[2][1]=0;A[2][2]=2;A[2][3]=1;A[2][4]=0;
A[2][5]=3;A[2][6]=2;A[2][7]=1;A[2][8]=0;A[2][9]=0;
A[2][10]=-1;A[2][11]=0;
A[3][0]=100;A[3][1]=1;A[3][2]=0;A[3][3]=1;A[3][4]=3;
A[3][5]=0;A[3][6]=2;A[3][7]=3;A[3][8]=5;A[3][9]=0;
A[3][10]=0;A[3][11]=-1;
break;
case3:
A[0][1]=5;A[0][2]=6;A[0][3]=23;A[0][4]=5;A[0][5]=24;
A[0][6]=6;A[0][7]=23;A[0][8]=5;
A[1][0]=100;A[1][1]=2;A[1][2]=1;A[1][3]=1;A[1][4]=1;
A[1][5]=0;A[1][6]=0;A[1][7]=0;
A[1][8]=0;A[1][9]=-1;
A[2][0]=150;A[2][1]=0;A[2][2]=2;A[2][3]=1;A[2][4]=0;
A[2][5]=3;A[2][6]=2;A[2][7]=1;A[2][8]=0;A[2][9]=0;
A[2][10]=1;A[2][11]=0;
A[3][0]=100;A[3][1]=1;A[3][2]=0;A[3][3]=1;A[3][4]=3;
A[3][5]=0;A[3][6]=2;A[3][7]=3;A[3][8]=5;A[3][9]=0;
A[3][10]=0;A[3][11]=-1;
break;
default:
cout<<"题号输入错误";
break;
}
}
intmain（）
{
intm,n,num,rk,rv,l,p;
inti,j,k,L[15];
doublev_min,mid,b,min,sum;
staticdoublex[10];
cout<<"请输入结构变量的个数:
"<cin>>m;
cout<<"请输入约束条件的个数:
"<cin>>n;
cout<<"请选择目标函数的极值（极大:
-1,极小：
1）:
"<cin>>p;
cout<<"请输入需求解问题的题号（1,2,3）:
"<cin>>num;
rk=0;
init_Array（num）;
for（i=0;i<=m;i++）
A[0][i]=A[0][i]*p;
for（i=1;i<=n;i++）
{
if（A[i][m+i]==1）
L[i]=m+i;
else
{
L[i]=m+n+i;
A[0][m+n+i]=M;
A[i][m+n+i]=1;
rk++;
}
}
if（rk!
=0）
{
for（j=1;j<=m+n+rk;j++）
A[n+1][j]=A[0][j];
for（j=1;j<=m+n+rk;j++）
{
sum=0;
for（i=1;i<=n;i++）
{
sum=sum+A[n+1][L[i]]*A[i][j];}
A[0][j]=A[n+1][j]-sum;//重求目标函数的系数
}
}
while
（1）
{
k++;
v_min=A[0][1];
l=1;
for（j=2;j<=m+n;j++）
{
if（A[0][j]{
v_min=A[0][j];
l=j;
}
}//求进基变量
rv=0;
for（i=1;i<=n;i++）
{
if（L[i]>m+n）
{
if（A[i][0]!
=0）
rv=1;
break;
}
}
if（v_min>=0）
{
if（rv==0）
{
for（i=1;i<=n;i++）
x[L[i]]=A[i][0];
min=function（x[1],x[2],x[3],x[4],x[5],x[6],x[7],x[8]）;
}//若目标函数中各变量的系数均为正，且该顶点中的人为变量等于0，则该顶点为最优解
else
{
cout<cout<<"noanswer"<return;
}
break;//若目标函数中各变量的系数均为正，且该顶点中的人为变量不等于0，则无解
}
else
{
b=A[1][l];
for（j=2;j<=n;j++）
{
mid=A[j][l];
if（mid>b）
b=mid;
}
if（b<0）
{
cout<cout<<"noanswer"<return;
}
for（i=1;i<=n;i++）
if（（A[i][0]/A[i][l]）>0）
{
b=A[i][0]/A[i][l];
k=i;
break;
}//必须为正才可以
for（i=1;i<=n;i++）
{
mid=A[i][0]/A[i][l];
if（（mid0））
{
b=mid;
k=i;
}//确定出基变量
}
L[k]=l;
for（i=0;i<=n;i++）
for（j=0;j<=m+n+rk;j++）
A[i+n+1][j]=A[i][j];
for（i=1;i<=n;i++）
{
if（i==k）
for（j=0;j<=n+m+rk;j++）
A[k][j]=A[k+n+1][j]/A[k+n+1][l];
else
{
for（j=0;j<=n+m+rk;j++）
A[i][j]=A[i+n+1][j]-A[k+n+1][j]*A[i+n+1][l]/A[k+n+1][l];}
}//进行线性变换
for（j=1;j<=m+n+rk;j++）
A[n+1][j]=A[0][j];
for（j=1;j<=m+n+rk;j++）
{
sum=0;
for（i=1;i<=n;i++）
{
sum=sum+A[n+1][L[i]]*A[i][j];
}
A[0][j]=A[n+1][j]-sum;
}//重新求线性方程组的系数
}
}
cout<<"********************************"<for（i=1;i<=8;i++）
{
cout<<"按截法"<x["<}
cout<<"总的边料最少为:
ObjectiveFunction="<cout<<"********************************"<return0;
}
运行结果如下：
（a）题：
（b）题：
（c）题：
动态规划法
第三题采用的动态规划法，动态规划从本质上讲是一种非线性规划方法，其核心是贝尔曼最优性原理，这些原理可以归结为一个基本的递推关系式，从而使决策过程连续的转移，并将一个多步最优控制问题化为多个一步最优控制问题，从而简化求解过程，动态规划的基本模型如下：
1.确定问题的决策对象。
2.对决策过程划分阶段。
3.对各阶段确定状态变量。
4.根据状态变量确定费用函数和目标函数。
5.建立各阶段状态变量的转移过程，确定状态转移方程。
多极决策过程的最优策略具有这样的性质，不论初始状态和初始决策如何，当把其中的任何一级和状态再作为初始级和初始状态时，其余决策对此必定也是一个最优决策。
3要求从城市1到12建立一条客运线。
各段路所能获得的利润如下表，注意：
从一个给定的城市出发只能直接到达某些城市。
问从城市1到12应走怎样的路线，才能获得最大利润。
到达城市
2
3
4
5
6
7
8
9
10
11
12
利润
出发城市
1
5
4
2
2
8
10
5
7
3
6
3
8
10
4
8
9
6
4
5
8
4
3
6
5
2
7
7
4
10
6
8
12
5
2
9
7
10
3
11
6
程序清单如下：
#include
#include
#include
usingnamespacestd;
constintN=11;
constintM=12;
staticintA[N][M];
intmain（）
{
A[1][2]=5;A[1][3]=4;A[1][4]=2;
A[2][5]=8;A[2][6]=10;A[2][7]=5;A[2][8]=7;
A[3][5]=6;A[3][6]=3;A[3][7]=8;A[3][8]=10;
A[4][5]=8;A[4][6]=9;A[4][7]=6;A[4][8]=4;
A[5][9]=8;A[5][10]=4;A[5][11]=3;
A[6][9]=5;A[6][10]=2;A[6][11]=7;
A[7][9]=4;A[7][10]=10;A[7][11]=6;
A[8][9]=12;A[8][10]=5;A[8][11]=2;
A[9][12]=7;
A[10][12]=3;
A[11][12]=6;//出发城市到到达城市的利润值
inttotal_0[3],total_1[4];
mapm;//城市号作为关联容器的键，下一站的的城市号为关联容器的值
vectordistant;//最优路径存放的向量
inti=0,j=0;
for（i=2;i<=4;i++）//第1站到第2站的利润值
{
total_0[i-2]=A[1][i];
m.insert（map:
:
value_type（i,1））;
}
for（i=5;i<=8;i++）//第2站到第3站的利润值
{
intmax=0,k;
for（j=2;j<=4;j++）
{
if（total_0[j-2]+A[j][i]>max）
{
max=total_0[j-2]+A[j][i];
k=j;
}
}
total_1[i-5]=max;
m.insert（map:
:
value_type（i,k））;
}
for（i=9;i<=11;i++）//第3站到第4站的利润值
{
intmax=0,k;
for（j=5;j<=8;j++）
{
if（total_1[j-5]+A[j][i]>ma
升级会员 