C++实现基因遗传驱动的列生成算法

明湖小虾

已于 2022-09-25 17:09:43 修改

阅读量437

点赞数

文章标签： c++ 算法 java

于 2022-09-24 18:12:37 首次发布

本文链接：https://blog.csdn.net/weixin_52466509/article/details/127028732

版权

列生成算法可适用于可用列进行线性化表示的非线性问题，其核心思想为使用对偶原理寻找最有可能改善原问题的解的列，最好的情况下，原问题为非线性的但对偶模型为线性模型，因此总能求解出一个好列，例如cutstock问题（卷纸切割问题）。

那么当对偶模型也为非线性问题无法使用线性规划的方法来求解新方案时，启发或元启发算法就可以派上用场了，可以求解出局部最优的进基列。

本文要求解的问题是指派、调度问题的一种，将m个工件分配到n个机床上进行加工，并且应当具有一定的顺序，使得最终的总用时最小的类似这种问题。

class MS {
public:
	int m, n;
	vector<vector<double>>p;//p_ij加工时间
	vector<double>D;//D_i
	vector<double>w;//w_j工作j的权重
public:
	void read();
	vector<oneG> getNewSbyGenetic(IloNumArray PI, IloNumArray RHO,int i);
	vector<bool> getNewSbyDP(IloNumArray PI, IloNumArray RHO);
	double getC(vector<bool> S,int i);
	void cplex();
	void cplexSolveMILP();
};

struct oneG {
public:
	double fitness;
	struct vector<bool>detail;
};

class Genetic {
public:
	int m, n;//机器数，作业数。
	int i;//子问题
	int number;//种群规模；
	vector<oneG>P;//现有种群。
	MS* ms;
	IloNumArray PI;
	IloNumArray RHO;
public:
	Genetic(MS* ms,int m, int n,int i,IloNumArray PI, IloNumArray RHO, int number=100) {
		this->m = m;
		this->n = n;
		this->number = number;
		this->ms = ms;
		this->i = i;
		this->PI = PI;
		this->RHO=RHO;
	}
	vector<oneG> initP();//种群初始化
	oneG randOne();//随机产生子代
	void cross();//交叉遗传；
};

cpp文件：

#include"MS.h"

bool sortFitness(oneG a, oneG b) { return (a.fitness < b.fitness); }
void MS::read() {//文件输入
	string tempString;
	ifstream readFile("data.txt");

	while (readFile >> tempString) {
		if (tempString == "机器数目:") {
			readFile >> m;//m个机器
		}
		if (tempString == "工作数目:") {
			readFile >> n;//n个工作
		}
		if (tempString == "加工时间:") {
			for (int i = 0; i < m; i++) {
				vector<double>pi;
				for (int j = 0; j < n; j++) {
					double pij;
					readFile >> pij;
					pi.push_back(pij);
				}
				p.push_back(pi);
			}
		}
		if(tempString=="中断时间:"){
			for (int i = 0; i < m; i++) {
				double di;
				readFile >> di;
				D.push_back(di);
			}
		}
		if (tempString == "工作权重:") {
			for (int j = 0; j < n; j++) {
				double wj;
				readFile >> wj;
				w.push_back(wj);
			}
		}
	}
	//std::cout << "数据输入完成" << endl;
}
void MS::cplexSolveMILP() {
	IloEnv env2;
	IloModel MILP(env2);
	IloArray<IloArray<IloNumVarArray>>X(env2);
	IloNumVarArray C(env2, n);
	for (int j = 0; j < n; j++) {
		string name = "C_" + to_string(j + 1);
		C[j] = IloNumVar(env2, 0, INFINITY,name.c_str());
	}
	IloRangeArray con1(env2);//第一组约束。
	IloArray<IloRangeArray> con2(env2);
	IloArray<IloRangeArray> con3(env2);
	IloArray<IloArray<IloRangeArray>> con4(env2);
	for (int j = 0; j < n; j++) {
		con1.add(IloRange(env2, 1, 1));
	}
	for (int r = 0; r < n; r++) {
		con2.add(IloRangeArray(env2));
		for (int i = 0; i < m; i++) {
			con2[r].add(IloRange(env2, -INFINITY, 1));
		}
	}
	for (int r = 0; r < n-1; r++) {
		con3.add(IloRangeArray(env2));
		for (int i = 0; i < m; i++) {
			con3[r].add(IloRange(env2, -INFINITY, 0));
		}
	}
	for (int i = 0; i < m; i++) {
		con4.add(IloArray<IloRangeArray>(env2));
		for (int j = 0; j < n; j++) {
			con4[i].add(IloRangeArray(env2));
			for (int r = 0; r < n; r++) {
				con4[i][j].add(IloRange(env2, -INFINITY, 0));
			}
		}
	}
	for (int i = 0; i < m; i++) {
		IloArray<IloNumVarArray>X_i(env2);
		for (int j = 0; j < n; j++) {
			IloNumVarArray X_ij(env2);
			for (int r = 0; r < n; r++) {
				string name = "X" + to_string(i + 1) + to_string(j + 1) + to_string(r + 1);
				X_ij.add(IloNumVar(env2, 0, 1, ILOBOOL,name.c_str()));
			}
			X_i.add(X_ij);
		}
		X.add(X_i);
	}
	for (int j = 0; j < n; j++) {
		IloExpr expr(env2);
		for (int i = 0; i < m; i++) {
			for (int r = 0; r < n; r++) {
				expr += X[i][j][r];
			}
		}
		con1[j].setExpr(expr);
	}

	for (int r = 0; r < n; r++) {
		for (int i = 0; i < m; i++) {
			IloExpr expr(env2);
			for (int j = 1; j < n; j++) {
				expr += X[i][j][r];
			}
			con2[r][i].setExpr(expr);
		}
		MILP.add(con2[r]);
	}
	for (int r = 0; r < n-1; r++) {
		for (int i = 0; i < m; i++) {
			IloExpr expr(env2);
			for (int j = 1; j < n; j++) {
				expr += X[i][j][r]-X[i][j][r+1];
			}
			con3[r][i].setExpr(expr);
		}
		MILP.add(con3[r]);
	}
	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) {
			for (int r = 0; r < n; r++) {
				IloExpr expr(env2);
				expr += p[i][j] * X[i][j][r];
				for (int l = 0; l < r - 1; l++) {
					for (int u = 1; u < n; u++) {
						expr += X[i][u][l] * p[i][u];
					}
				}
				for (int v = 0; v < r; v++) {
					for (int q = r + 1; q < n; q++) {
						for (int k = 0; k < n; k++) {
							expr += D[i] * pow(1 - D[i], v - 1) * p[i][k] * X[i][k][q];
						}
					}
				}
				expr -= C[j]+M*(1-X[i][j][r]);
				con4[i][j][r].setExpr(expr);
			}
			MILP.add(con4[i][j]);
		}
	}
	IloObjective obj;
	IloExpr expr2(env2);
	for (int j = 0; j < n; j++) {
		expr2 += w[j]*C[j];
	}
	obj = IloMinimize(env2, expr2);
	MILP.add(obj);
	MILP.add(con1);
	
	IloCplex MILPSolver(MILP);
	MILPSolver.exportModel("MILPmodel.lp");
	MILPSolver.solve();
	cout<<"求解结果:"<< MILPSolver.getStatus()<<endl;
	cout << "最优解:" << MILPSolver.getObjValue()<<endl;
	IloNumArray val(env2);
	MILPSolver.getValues(C,val);
	cout << val<<endl;
	for (int i = 0; i < m; i++) {
		cout << "机器" + to_string(i + 1) + "上的工作:";
		for (int j = 0; j < n; j++) {
			IloNumArray val(env2);
			MILPSolver.getValues(X[i][j], val);
			for (int r = 0; r < n; r++) {
				if (val[r] > 0.5) {
					cout << j << "工作在"<<r<<"位置、";
				}
			}
		}
		cout << endl;
	}
	/*for (int i = 0; i < m; i++) {
		cout << "机器" + to_string(i + 1) + "上的工作:";
		for (int j = 0; j < n; j++) {
			IloNumArray val(env2);
			MILPSolver.getValues(X[i][j], val);
			cout << val<<endl;
		}
	}*/
}
void MS::cplex() {
	IloEnv env;
	try {
		//主问题
		IloModel msOpt(env);
		IloObjective rollsUsed = IloAdd(msOpt, IloMinimize(env));
		IloRangeArray constraints1(env);//第一组约束。
		IloRangeArray constraints2(env);
		for (int j = 0; j < n; j++) {
			constraints1.add(IloRange(env, 1, 1));
		}
		for (int i = 0; i < m; i++) {
			constraints2.add(IloRange(env, -INFINITY, 1));
		}
		//env.out() << constraints1.getSize();
		IloArray<IloNumVarArray> S(env);
		vector<vector<bool>>initSCs{ {0,0,0,1,0,1,0,1,0,0},{1,1,1,0,0,0,0,0,1,1},{0,0,0,0,1,0,1,0,0,0} };
		for (int i = 0; i < m; i++) {
			S.add(IloNumVarArray(env));
			//初始化列
			IloNumArray initS(env, n);
			//
			vector<bool>initSC=initSCs[i];
			//
			for (int ii = 0; ii < n; ii++) {
				initS[ii] = initSC[ii];
			}
			/*精确模型最优解
			机器1上的工作:3、5、7、
			机器2上的工作 : 0、1、2、8、9、
			机器3上的工作 : 4、6、*/
			double sj=getC(initSC,i);
			IloExpr expr(env);
			for (int ii = 0; ii < 1; ii++) {
				string name = "lambda" + to_string(i+1) + to_string(ii+1);
				S[i].add(IloNumVar(rollsUsed(sj) + constraints1(initS)));
				//cout << S[i]<<endl;
				S[i][ii].setName(name.c_str());
				expr += S[i][ii];
			}
			constraints2[i].setExpr(expr);
		}
		
		msOpt.add(constraints1);
		msOpt.add(constraints2);
		IloCplex msSolver(msOpt);
		//msSolver.exportModel("model.lp");
		bool Generate = 1;
		while (Generate) {
			Generate = 0;
			msSolver.exportModel("model.lp");
			msSolver.solve();
			IloNumArray price1(env, n);
			IloNumArray price2(env, m);
			IloNumArray newSolutions(env, n);
			//输出对偶值。
			cout << msSolver.getStatus() << endl;
			cout << msSolver.getObjValue() << endl;;
			for (int j = 0; j < n; j++) {
				price1[j] = msSolver.getDual(constraints1[j]);
				cout << price1[j] <<" ";
			}
			cout << endl;
			for (int i = 0; i < m; i++) {
				price2[i] = msSolver.getDual(constraints2[i]);
				cout << price2[i] << " ";
			}
			cout << endl;
			
			//子问题
			for (int i = 0; i < m; i++) {
				vector<oneG> Gs = getNewSbyGenetic(price1, price2, i);
				for (int l = 0; l < Gs.size(); l++) {
					oneG G = Gs[l];
					vector<bool> newSC = G.detail;
					for (int j = 0; j < newSC.size(); j++) {
						newSolutions[j] = newSC[j];

					}
					cout << newSolutions;
					double sj = getC(newSC, i);
					cout << G.fitness << endl;
					if (G.fitness < -RC_EPS) { Generate = 1; };
					IloNumArray newCons2(env, m);
					for (int iii = 0; iii < m; iii++) {
						newCons2[iii] = 1;
					}
					string name = "lambda" + to_string(i + 1) + to_string(S[i].getSize() + 1);
					S[i].add(IloNumVar(rollsUsed(sj) + constraints1(newSolutions) + constraints2(newCons2)));
					S[i][S[i].getSize() - 1].setName(name.c_str());
				}
			}
		}
		for (int i = 0; i < S.getSize(); i++) {
			msOpt.add(IloConversion(env, S[i], ILOBOOL));
		}
		msSolver.solve();
		msSolver.exportModel("model.lp");
		cout << "Solution status: " << msSolver.getStatus() << endl;
		cout << "最优解:" << msSolver.getObjValue() << endl;
		for (int i = 0; i < S.getSize(); i++) {
			IloNumArray val(env);
			msSolver.getValues(S[i], val);
			cout << val<<endl;

		}
	}
	catch (IloException& ex) {
		cerr << "Error: " << ex << endl;
	}
	catch (...) {
		cerr << "Error" << endl;
	}

	env.end();

}
vector<oneG> MS::getNewSbyGenetic(IloNumArray PI, IloNumArray RHO,int i) {
	vector<oneG> solutions;
	Genetic genetic(this,m, n,i, PI, RHO);
	//cout <<"第"<<i<< "问题初始化种群" << endl;
	genetic.initP();
	//cout << "开始遗传算法" << endl;
	for (int k = 0; k < 50; k++) {
		genetic.jiaocha();
	}
	for (int k = 1;k < genetic.P.size()+1; k++) {
		if (genetic.P[k].fitness != genetic.P[k - 1].fitness) {
			solutions.push_back(genetic.P[k - 1]);//添加新方案组
		}

		if (solutions.size() > 10|| genetic.P[k].fitness>0)break;
	}
	//cout << "找到可行新方案" << endl;
	return solutions;
}
vector<bool> MS::getNewSbyDP(IloNumArray PI, IloNumArray RHO) {
	vector<bool>solution;
	return solution;
}
double MS:: getC(vector<bool> S,int i) {
	//计算方案的代价
	double csi = 0;
	double t=0;
	vector<int>work;
	for (int j = 0; j < n; j++) {
		if (S[j] != 0) {
			work.push_back(j);
		}
	}
	for (int j = 0; j < work.size(); j++) {
		t += p[i][work[j]] * pow(1 - D[i], i);
		for (int jj = j+1; jj < work.size(); jj++) {
			t += p[i][work[jj]] * D[i]*pow(1-D[i],i);
		}
		csi += t * w[work[j]];
	}
	return csi;
}
vector<oneG> Genetic::initP(){
	//初始化种群
	vector<oneG>Q;
	for (int k = 0; k < number; k++) {
		oneG np = randOne();
		Q.push_back(np);
	}
	P = Q;
	return P;
}
oneG Genetic::randOne() {
	oneG randO;
	for (int k = 0; k < n; k++) {
		randO.detail.push_back(0);
	}
	for (int k = 0; k < n; k++) {
		int temp= (rand() % (n + 1));
		double prob = double(rand() % 10)/10;
		if (prob < 0.3) {
			randO.detail[temp]=1;
		}
	}
	//计算适应度
	randO.fitness= ms->getC(randO.detail, i);
	//根据对偶值计算
	for (int j = 0; j < n; j++) {
		randO.fitness -=randO.detail[j]*PI[j];
	}
	for (int j = 0; j < m; j++) {
		randO.fitness -= RHO[j];
	}
	return randO;
}
void Genetic::cross() {
	double sum = 0;
	vector<oneG>Q;//新种群
	//cout << "交叉遗传中" << endl;
	for (int k = 0; k < P.size(); k++) {
		for (int ii = k+1; ii < P.size(); ii++) {
			oneG ne;
			int point = (rand() % (n + 1));
			for (int j = 0; j < n; j++) {
				if (j < point) {
					ne.detail.push_back(P[k].detail[j]);
				}
				else {
					ne.detail.push_back(P[ii].detail[j]);
				}
			}
			ne.fitness = ms->getC(ne.detail, i);
			for (int j = 0; j < n; j++) {
				ne.fitness -= ne.detail[j] * PI[j];
			}
			for (int j = 0; j < m; j++) {
				ne.fitness -= RHO[j];
			}
			Q.push_back(ne);
		}
	}
	Q.insert(Q.end(), P.begin(), P.end());//加父本
	sort(Q.begin(), Q.end(), sortFitness);//排序
	P.resize(0);
	for (int k = 1; k < Q.size()+1; k++) {//去重复种群
		if (Q[k].fitness != Q[k - 1].fitness) {
			P.push_back(Q[k - 1]);
		}
		if (P.size() > number)break;
	}
	//cout << "交叉完成" << endl;
}

主程序：

#include"MS.h"

ILOSTLBEGIN

using namespace std;

void main() {
	MS ms;
	ms.read();
	int a = 0;
	cout << "请选择求解方式：1.MILP精确求解；2.基因遗传列生成求解。" << endl;
	cin >> a;
	if (a == 2) ms.cplex();
	else ms.cplexSolveMILP();
}

明湖小虾

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
C++实现基因遗传驱动的列生成算法

列生成算法可适用于可用列进行线性化表示的非线性问题，其核心思想为使用对偶原理寻找最有可能改善原问题的解的列，最好的情况下，原问题为非线性的但对偶模型为线性模型，因此总能求解出一个好列，例如cutstock问题（卷纸切割问题）。那么当对偶模型也为非线性问题无法使用线性规划的方法来求解新方案时，启发或元启发算法就可以派上用场了，可以求解出局部最优的进基列。本文要求解的问题是指派、调度问题的一种，将m个工件分配到n个机床上进行加工，并且应当具有一定的顺序，使得最终的总用时最小。
复制链接

扫一扫