割平面法编程
云南大学数学与统计学实验教学中心
实验报告
一、实验目的
编程实现割平面整数线性规划问题，加深对割平面法以及单纯形法的理解和掌握
二、实验环境
VS2010（C++）
三、实验内容
在对偶单纯形的基础上，考虑x 的不可分性，也就是不能为小数。这样就需要我们在求的最优解得基础上，保证x的取值为整数。
经过学习，有两种方法可以用来求解：分支定界法，以及割平面法。
割平面法就是通过一定的约束，将非整数解割去，在其凸边上得到相应的整数解。
首先不考虑变量xi是整数这一条件，但是增加线性约束条件（称为割平面），使得原来可行解区域切去一部分，但没有切去任何整数可行解。该方法是指出怎样找到适当的割平面，使得切割后得到可行解区域的某顶点（称极点，对应基本可行解）恰好是整数线性规划问题的最优解
四、实验过程
A、割平面法的算法思想（课本割平面法部分由具体算法P121）
1.给定整数线性规划（ILP）的实例，用原始单纯形算法解决ILP的松弛线性规划（LP），
得到LP的最优基本可行解x，其基为B，最后单纯形表中的第i个方程为（i=1,2,…,m)
xi + ∑k aik xk = bi (5.1) 由于方程5.1中的变量是非负的，则
∑k [aik]xk ≤∑k aik xk (5.2)
从而方程(5.1)可变为
xi + ∑k [aik]xk ≤bi （5.3 因为在ILP中x的限制为整数，从而(5.3)的左边为整数，进而(5.3)的右边也为整数，
于是得到
xi + ∑k [aik]xk ≤[bi] (5.4) 用(5.1)减去(5.4)，得到
∑k（aik - [aik]）xk ≥bi - [bi] (5.5)
令fik=aik - [aik]和fi=bi - [bi]，这里fik（fi）称为aik
（bi）小数部分，满足
0 ≤fik< 1 和0 ≤fi < 1 (5.6)
从而由(5.5)和 (5.6)得到下面的式子然后再按照原单纯形算法经行求解。重复上述步骤。
2. 直到得到的解为整数。才停止运算。
B 、编译运行程序，输入约束方程的系数矩阵A ，常矩阵B,价值系数C,得到最优解
为了保证得到的算法能求得解，并且求的的最优解是对的。这里选取了P118页的例题3。
1、 约束条件
12
12121212max 134,0,z x x x x x x x x x x =+-+≤??+≤??
≥???为整数
2、 整理化简后，得到初始表如下：
Cj
1 1 0 0 CB
XB b x1 x2 x3 x4 0 x3 1 -1 1 1 0 0
x4
4
3
1
1
D 、运行结果如下图：
1）代数上述表中的值：
a 、初始化后的结果：这时候，它是按照对偶单纯形来求解的
b 、对偶单纯形优化后的结果：
得到最优解，但是此时得到的解为小数而非整数解
c、运用割平面法经行优化改进
割平面法添加条件后出事后以及第一次优化的结果
割平面法的第二次第三次优化：
d、最后得到的结果：
最后得到的结果与书上例题所有的结果完全吻合。所以所编写的算法为割平面法。
五、实验总结
通过对割平面算法编程进一步加强了对本章知识的掌握，当然这不局限于割平面法，对分支定界法也有进一步的了解。
六、源代码
这里需要说明的是。这个算法是在对偶单纯形算法的基础上修来得来的。
源代码如下：
#include
#include
using namespace std;
#define M -10000
class DanChun{
public:
double A[100][120],B[30]; //创建一个二维数组，用来存储xi值以及b的值
double A1[100][120],B1[30],C1[110]; //创建一个二维数组，用来存储xi值以及b的值
double C[110],CB[40];//定义一个C矩阵用来存储x 的价值系数
double CZ[110];//CZ用来存储cj-zj的值
double sita[30];//用来存储sita的值
int XB[30];//用来存XB的值
double X[100];//用来存储最后计算的X的值
double w[100]; //Cj-Zj/X[i]的替换。
public:
int n,m;
double favl;
double XZ[100]; // 得到最优解以及xi 的值
double RA[100][120]; // 得到最优解时A的系数
void init();//1.将不等式约束的各矩阵进行赋值
bool FindE(int i, int j);//2.寻找A矩阵中的基向量
void FindCB(); //3.寻找有用价值系数，Cj-Zj需要用到的
int Cj_Zj(); //4.1求出cj-zj的值
int C_sita(int k);//4.2计算sita的值
void Com_DC();//4.用来计算单纯形法的主函数
void Com_E(int l, int k);//4.3求其单位向量
bool P_CjZj();//4.4
int Bbool(double b[], int k); //判断b中是否有负值存在int TH(int l); //b有负值是，进行替换操作
void Display();
void Display1();
void Setinit(int m, int n, double a[][120], double *b);
void Reback();
/*----*/
void Buffer(double b,int k);
void start();
void Gepingmianfa();
bool panduan(double x[100]);
};
//增加条件后，传过来的行数，列数，以及初始值
//a为约束条件系数
//c 为价值系数
//b 为约束系数
void DanChun::init(){
int i=0,j=0;
cout<<"请输入 m,n ：";
cin>>m;
cin>>n;
cout<<"请输入线性约束矩阵："< for(i=0; i for(j=0;j if(j cin>>A[i][j]; } else { A[i][j]=0; if(j==n+i){ A[i][j]=1; } } A1[i][j] = A[i][j]; } } cout< cout<<"请输入B的值："< for(i=0;i cin>>B[i]; B1[i] = B[i]; CB[i]=0; } cout< cout<<"请输入价值系数 C："< for(i=0;i if(i cin>>C[i]; } else{ C[i] = M; } X[i]=0; C1[i] = C[i]; } n = m+n; } void DanChun::Display(){ int i,j; cout<<"价值系数 C："; for(i=0;i cout< } cout< cout<<"B的值："; for(i=0;i cout< } cout< cout<<"线性约束矩阵A："< for(i=0; i for(j=0;j printf("%0.3lf\t",A[i][j]); } cout< } cout< cout<<"价值系数 CB："; for(i=0;i cout< } cout< } bool DanChun::P_CjZj(){ int j; for(j=0;j if(CZ[j]>0){ return true; } } return false; } int DanChun::TH(int l){ int i=0, k=0; double min = 1000; for(i=0; i w[i] = 0; } for(i=0; i if(A[l][i]>=0){ w[i] = -M; }else{ w[i] = CZ[i]/A[l][i]; } } for(i=0; i if(min>w[i]){ min = w[i]; k = i; } } return k; } int DanChun::Bbool(double b[], int k){ int i=0, l=0, s=0; double min = 0; for(i=0; i if(min>b[i]){ min = b[i]; l = i; //找出行下标 } } if(min< 0){ return l; } return -1; } void DanChun::Display1(){ int i,j=0,p=0; double favl=0; for( i=0;i X[XB[i]]=B[i]; } for(i=0;i favl += C[i]*X[i]; } cout<<"当前为最优解：favl = "< cout<<"有最优解时，x的取值如下："< for(i=0;i j=i; cout<<"X["< } } void DanChun::Com_DC(){ int i,j=0,k,l,p=0; bool b; double favl=0; do{ b =true; k=Cj_Zj();//得到cj-zj的大于0的max if(k == -1){ if( (l=Bbool(B,m)) >=0){ k = TH(l); b = false; } else { Display1(); return; } } if(b){ l=C_sita(k); if(l==-1){ cout<<"无解！"< return; } } CB[l]=C[k]; XB[l]=k; Com_E(l,k); p++; cout<<"----------------------------第"< cout< for(i=0;i cout< } cout< Display(); cout<<"----------------------------第"< }while(P_CjZj() || !b); } void DanChun::Com_E(int l, int k){ int i,j; B[l]=B[l]/A[l][k]; for(j=0;j if(j==k) continue; A[l][j]= A[l][j]/A[l][k]; } A[l][k]=1; for(i=0;i if(i==l) continue; for(j=0;j if(j==k) continue; A[i][j]=A[i][j]-A[i][k]*A[l][j]; } B[i]=B[i]-A[i][k]*B[l]; A[i][k]=0; } } int DanChun::Cj_Zj(){ int i,j,k=0; double sum=0,max=0; for(j=0;j for(i=0;i sum+=CB[i]*A[i][j]; } CZ[j]=C[j]-sum; sum=0; } for(j=0;j if(CZ[j]>max && CZ[j]>0){ max=CZ[j]; k=j; } } if(max==0){ return -1; } return k; } int DanChun::C_sita(int k){ int i,j=0; double min=-M; for(i=0;i if(A[i][k]>0){ sita[i]=B[i]/A[i][k]; j++; }else{ sita[i]=-M; } } if(j==0){ //表示没有一个A[i][k]是大于0的 return -1; } for(i=0;i if(sita[i] min=sita[i]; k=i; } } return k; } bool DanChun::FindE(int i,int j){ //判断是否是单位列int k=0; while(k if(A[k][j] !=0 && k!=i){ return false; }else{ k++; if(k==i){ continue; } } } return true; } void DanChun::FindCB(){ int i,j; for(i=0; i for(j=0;j if(A[i][j]==1){ if(FindE(i,j)==true){ //如果是单位列，怎选择为其对应的价值系数 CB[i]=C[j]; XB[i]=j; break; } } } } } /*------------------------割平面法------------------------------------*/ void DanChun::Setinit(int m1, int n1, double a[][120],double b[]){ //将值传入int i,j; m = m1; n = n1; //约束系数 for(i=0; i for(j=0;j if(j A[i][j] = a[i][j]; } else { A[i][j]=0; if(j==n+i){ A[i][j]=1; } } } } for(i=0;i B[i] = b[i]; CB[i]=0; } for(i=0;i if(i C[i] = C1[i]; C[i] = M; } X[i]=0; } } void DanChun::Reback(){ int i,j; //XZ[100]; 得到当前最优解 //以及 xi 的值 XZ[0] = favl; for(i=0;i XZ[i+1]= X[i]; } //得到想x系数的矩阵A for(i=0; i for(j=0;j RA[i][j] = A[i][j]; } } } void DanChun::Buffer(double b,int k){ int i,j; double X[100]; Reback(); //对应的第k行向下取整 for(i=0; i X[i] = floor(RA[k][i]); } //将X[i]中的值加到原来的约束中，即增加一行for(i=0; i A1[m][i] = X[i] - RA[k][i]; } //相应的b的系数增加一个 B1[m] = b; //增加一行 int m1 = m + 1; //增加一列 int n1 = n + 1; double s = 0,t =0; for(i=0; i if(i==(n-m)){ s = C1[i]; C1[i] = 0; for(j=i+1; j< n1; j++){ C1[j] = s; t = C1[j]; s = t; } } for(j=n-1; j>=n-m;j--){ for(i=0; i A1[i][j+1] = A1[i][j]; } } for(i=0; i if(i A1[i][n-m] = 0; }else{ A1[i][n-m] = 1; } } Setinit(m1,n1,A1,B1); FindCB(); Display(); Com_DC(); } void DanChun::start(){ init(); FindCB(); Display(); Com_DC(); Gepingmianfa(); } void DanChun::Gepingmianfa(){ int i=0,j=0; bool b; double XZ1[100]; //得到XZ中的值（max z ,xi） double eps = 0.000001; double Fb; //向下取整floor得到的值 for(i = 0; i XZ1[i] = XZ[i]; } i = 0; int q = m; while(i //if(Xlocation[i] - floor(Xlocation[i]) > eps){ if(X[i] - floor(X[i]) > eps){ Fb = floor(X[i]); Fb = Fb - X[i]; Buffer(Fb,i); b = panduan(X); }else{ b=false; if(i==q-1){ b = true; } } if(b){ cout<<"此时有最优解！"< break; } i++; } } bool DanChun::panduan(double x[]){ int i; bool b; double eps = 0.000001; for(i=0;i if(x[i] - floor(x[i]) > eps && (ceil(x[i]) - x[i]) > eps ){ b = false; return b; }else { b = true; } } return b; } void main(){ DanChun D; D.start(); }