凸二次规划解法（旋转法）——补充代码

根据之前的两篇博客，理解了凸二次规划解法（旋转法）的计算过程和步骤后，根据自己手算的计算过程和步骤，用机器表示。

C语言代码的步骤和课文上的步骤是一样：

C语言凸二次规划的计算结果如下：（例子为上篇博客的例子）上篇博客：点击打开链接

c语言代码如下：

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "vector"
using namespace std;


#define MAXITER 1000  //最大迭代次数
#define INF 99999
#define varNum 2  //变量数
#define resNum1 0  //约束条件：非基等式
#define resNum2 2 //约束条件：基等式
#define COEF  1   //规划的系数


typedef vector<double> dim1Vector;   //1维浮点容器
typedef vector<dim1Vector> dim2Vector;  //2维浮点容器


void cal_LP(dim2Vector E, dim1Vector DEV, dim1Vector X);    //解线性规划
dim2Vector Cal_Spin(int laber_i, int laber_j, dim2Vector E, dim1Vector DEV);  //旋转运算


void main()
{
int i, j;
dim1Vector tempDim1;
double tempNum;
dim2Vector H;   //正定矩阵H
dim1Vector C;   
dim2Vector A1;   //约束条件：非基等式
dim2Vector A2;   //约束条件：基等式
dim1Vector B;    //约束条件结果
dim1Vector X;    //X的约束条件
dim1Vector tempE(varNum+resNum1+resNum2, 0);
dim2Vector E(varNum+resNum1+resNum2, tempE);    //初始基向量矩阵
dim1Vector DEV;  //偏差
double sum;

printf("输入 %d阶半正定矩阵H：\n", varNum);
for(i=0; i<varNum; i++)
{
tempDim1.clear();
for(j=0; j<varNum; j++)
{
scanf("%lf", &tempNum);
tempDim1.push_back(tempNum*COEF);
}

H.push_back(tempDim1);
}


printf("输入矩阵C：\n");
for(i=0; i<varNum; i++)
{
scanf("%lf", &tempNum);
C.push_back(tempNum);
}


printf("约束条件：非基等式A1\n");
for(i=0; i<resNum1; i++)
{
tempDim1.clear();
for(j=0; j<varNum; j++)
{
scanf("%lf", &tempNum);
tempDim1.push_back(tempNum);
}


A1.push_back(tempDim1);
}


printf("约束条件：基等式A2\n");
for(i=0; i<resNum2; i++)
{
tempDim1.clear();
for(j=0; j<varNum; j++)
{
scanf("%lf", &tempNum);
tempDim1.push_back(tempNum);
}

A2.push_back(tempDim1);
}




printf("约束条件结果B\n");
for(i=0; i<resNum1+resNum2; i++)
{
scanf("%lf", &tempNum);
B.push_back(tempNum);
}




printf("X的约束条件\n");
for(i=0; i<varNum+resNum1+resNum2; i++)
{
scanf("%lf", &tempNum);
X.push_back(tempNum);
}


//获取初始矩阵E
for(i=0; i<H.size(); i++)
for(j=0; j<H[i].size(); j++)
E[i][j] = H[i][j];


for(i=0; i<A1.size(); i++)
for(j=0; j<A1[i].size(); j++)
E[i+varNum][j] = A1[i][j];


for(i=0; i<A1.size(); i++)
for(j=0; j<A1[i].size(); j++)
E[j][i+varNum] = -A1[i][j];


for(i=0; i<A2.size(); i++)
for(j=0; j<A2[i].size(); j++)
E[i+varNum+resNum1][j] = A2[i][j];


for(i=0; i<A2.size(); i++)
for(j=0; j<A2[i].size(); j++)
E[j][i+varNum+resNum1] = -A2[i][j];




//获取偏差DEV
for(i=0; i<C.size(); i++)
DEV.push_back(C[i]);


for(i=0; i<B.size(); i++)
DEV.push_back(B[i]);

//初始化偏差
for(i=0; i<DEV.size(); i++)
{
sum = 0;
for(j=0; j<varNum; j++)
sum += E[i][j]*X[i];
DEV[i] = sum - DEV[i];
}


printf("初始情况：\n");

for(i=0; i<E.size(); i++)
{
for(j=0; j<E[i].size(); j++)
printf("%0.2lf  ", E[i][j]);
printf("%0.2lf\n", DEV[i]);
}


//开始解线性规划
cal_LP(E, DEV, X);
}




//解线性规划
void cal_LP(dim2Vector E, dim1Vector DEV, dim1Vector X)
{
int i, j, m, n;
int laber_i, laber_j;
dim2Vector calResult;
double maxNum;
double count;
dim1Vector Col;
dim1Vector Row(DEV.size(), INF);
double col2row;


for(i=0; i<E[0].size(); i++)
Col.push_back(i);

//非基等式转为基等式
for(i=0; i<resNum1; i++)
{
maxNum = 0;
count = 0;
//正数
if(DEV[i+varNum]>0)
{
for(j=0; j<E[i+varNum].size(); j++)
if(E[i+varNum][j]<0 && maxNum<fabs(E[i+varNum][j]))
{
maxNum = fabs(E[i+varNum][j]);
laber_i = i+varNum;
laber_j = j;
count++;
}

if(count==0)
{
printf("方程组无解!!\n");
exit(0);
}
}

//负数
if(DEV[i+varNum]<0)
{
for(j=0; j<E[i+varNum].size(); j++)
if(E[i+varNum][j]>0 && maxNum<fabs(E[i+varNum][j]))
{
maxNum = fabs(E[i+varNum][j]);
laber_i = i+varNum;
laber_j = j;
count++;
}

if(count==0)
{
printf("方程组无解!!\n");
exit(0);
}
}

//0
if(DEV[i+varNum]==0)
{
for(j=0; j<E[i+varNum].size(); j++)
if(E[i+varNum][j]!=0 && maxNum<fabs(E[i+varNum][j]))
{
maxNum = fabs(E[i+varNum][j]);
laber_i = i+varNum;
laber_j = j;
count++;
}

if(count==0)
{
printf("方程组无解!!\n");
exit(0);
}
}

col2row = Col[laber_j];
Col[laber_j] = Row[laber_i];
Row[laber_i] = col2row;

//计算旋转
calResult = Cal_Spin(laber_i, laber_j, E, DEV);

for(m=0; m<calResult.size()-1; m++)
for(n=0; n<calResult[m].size(); n++)
E[m][n] = calResult[m][n];

for(n=0; n<calResult[m].size(); n++)
  DEV[n] = calResult[m][n];

}



//开始迭代
int iter = 0;
bool flag = true;
double minNum;


while(iter<MAXITER)
{
minNum = 0;
count = 0;


for(i=0; i<DEV.size(); i++)
if(DEV[i]<0 && minNum>DEV[i])
{
minNum = DEV[i];
laber_i = i;
count++;
}


//迭代结束判断
if(count==0)  
flag = false;

if(!flag)
break;


//寻找旋转点
if(E[laber_i][laber_i]>0)
{
laber_j = laber_i;
count++;
}


else
{
for(j=0; j<E[laber_i].size(); j++)
if(E[laber_i][j]>0)
{
laber_j = j;
count++;
break;
}
}

if(count==0)
{
printf("方程组无解!!\n");
exit(0);
}


col2row = Col[laber_j];
Col[laber_j] = Row[laber_i];
Row[laber_i] = col2row;


//开始旋转
calResult = Cal_Spin(laber_i, laber_j, E, DEV);

for(m=0; m<calResult.size()-1; m++)
for(n=0; n<calResult[m].size(); n++)
E[m][n] = calResult[m][n];

for(n=0; n<calResult[m].size(); n++)
  DEV[n] = calResult[m][n];


iter++;
}


//计算X的解
for(i=0; i<Row.size(); i++)
{
if(Row[i]==INF)
continue;

X[Row[i]] = X[Row[i]]+DEV[i];
}


printf("\n迭代次数：%d\n", iter);

printf("旋转矩阵结果：\n");

for(i=0; i<E.size(); i++)
{
for(j=0; j<E[i].size(); j++)
printf("%0.2lf  ", E[i][j]);
printf("%lf\n", DEV[i]);
}


printf("线性规划的解为：\n");
for(i=0; i<varNum; i++)
printf("X%d = %lf  ", i+1, X[i]);
printf("\n");
}






//旋转运算
dim2Vector Cal_Spin(int laber_i, int laber_j, dim2Vector E, dim1Vector DEV)
{
int i, j;
dim1Vector temp(E[0].size(), 0);
dim2Vector dst(E.size()+1, temp);
dim2Vector E_TEMP;
dim1Vector DEV_TEMP(DEV.size(), 0);


E_TEMP = E;


//偏差旋转计算
DEV_TEMP[laber_i] = -DEV[laber_i]/E[laber_i][laber_j];
for(i=0; i<DEV.size(); i++)
{
if(i==laber_i)
continue;
DEV_TEMP[i] = DEV[i]-(E[i][laber_j]/E[laber_i][laber_j])*DEV[laber_i];
}




//矩阵偏差计算
dst[laber_i][laber_j] = INF;


for(j=0; j<E[laber_i].size(); j++)
{
if(dst[laber_i][j]==INF)
continue;
E_TEMP[laber_i][j] = -E[laber_i][j]/E[laber_i][laber_j];
dst[laber_i][j] = INF;
}


for(i=0; i<E.size(); i++)
{
if(dst[i][laber_j]==INF)
continue;
E_TEMP[i][laber_j] = E[i][laber_j]/E[laber_i][laber_j];


dst[i][laber_j] = INF;
}


for(i=0; i<E.size(); i++)
for(j=0; j<E[i].size(); j++)
{
if(dst[i][j]==INF)
continue;
E_TEMP[i][j] = E[i][j]-E[laber_i][j]*(E[i][laber_j]/E[laber_i][laber_j]);
}




E_TEMP[laber_i][laber_j] = 1/E[laber_i][laber_j];




//赋值、return
for(i=0; i<E_TEMP.size(); i++)
for(j=0; j<E_TEMP[i].size(); j++)
dst[i][j] = E_TEMP[i][j];


for(j=0; j<DEV_TEMP.size(); j++)
dst[i][j] = DEV_TEMP[j];

return dst;
}