根据之前的两篇博客,理解了凸二次规划解法(旋转法)的计算过程和步骤后,根据自己手算的计算过程和步骤,用机器表示。
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;
}