Raw File
new_ilp_solver.cpp
#include <bits/stdc++.h>
#include "gurobi_c++.h"
#include "Instance.h"

using namespace std;

struct node {
	int x, y, t, d;
	// (robot, xpos, ypos, time, dir)
	// dir == 5 is undecided
	node() {}
	node(int x, int y, int t, int d) : x(x), y(y), t(t), d(d) {}
	const bool operator<(const node& o) const {
		return make_tuple(t, x, y, d) < make_tuple(o.t, o.x, o.y, o.d);
	}
};

struct state {
	int r;
	node n;
	// (robot, xpos, ypos, time, dir)
	// dir == 5 is undecided
	state() {}
	state(int r, int x, int y, int t, int d) : r(r), n(x, y, t, d) {}
	const bool operator<(const state& o) const {
		if (r != o.r) return r < o.r;
		return n < o.n;
	}
};

int main() {
	string filename;
	cin >> filename;
	filename = remove_ext(filename);

	GRBEnv env(true);
	env.start();
	GRBModel model(env);

	map<state, GRBVar> Y;
	map<node, GRBLinExpr> Z;
	instance ins;
	ins.read(filename);


	const int minx = -5, maxx = 15, miny = -5, maxy = 15;
	auto in_bounds = [&](int x, int y) {
		return minx <= x && x <= maxx && miny <= y && y <= maxy;
	};

	// flow constraints
	int maxT = 30;
	
	for (int t = maxT; t >= 0; t--) {
		cerr << t << endl;
		for (int x = minx; x <= maxx; x++) {
			for (int y = miny; y <= maxy; y++) {
				// Each robot moves in the direction they've decided in
				for (int d = 0; d <= 5; d++) {
					vector<GRBVar> terms;
					vector<double> coeffs;
					for (int r = 0; r < ins.n; r++) {
						state s(r, x, y, t, d);
						// stringstream ss;
						// ss << "(" << x << "," << y << "," << t << "," << d << ")";
						Y[s] = model.addVar(0, 1, 0, GRB_BINARY);//, ss.str());
						
						if (d != 5 && t < maxT && in_bounds(x + dx[d], y + dy[d])) {
							state ns(r, x + dx[d], y + dy[d], t + 1, 5);
							// stringstream ss;
							// ss << "Adding constraint. Most move to " << x + dx[d] << " " << y + dy[d] << " " << d 
								 // << " from " << x << " " << y << " " << d;
							assert(Y.count(ns));
							model.addConstr(Y[s] <= Y[ns]);//, ss.str());

							// assert(Z.count(s.n));
							terms.push_back(Y[s]);
							coeffs.push_back(1);
						}
					}

					// Each state can only be occupied by at most 1 robot
					node curr(x, y, t, d);
					Z[curr].addTerms(&coeffs[0], &terms[0], terms.size());
					model.addConstr(Z[curr] <= 1);
				}
			}
		}
	}

	// Each square must take in a robot if its active
	for (int t = maxT; t >= 1; t--) {
		cerr << t << endl;
		for (int x = minx; x <= maxx; x++) {
			for (int y = miny; y <= maxy; y++) {
				// Each robot moves in the direction they've decided in
				for (int r = 0; r < ins.n; r++) {
					state s(r, x, y, t, 5);
					GRBLinExpr expr;
					for (int d = 0; d <= 4; d++) {
						if (in_bounds(x - dx[d], y - dy[d])) {
							state ns(r, x - dx[d], y - dy[d], t - 1, d);
							assert(Y.count(ns));
							expr += Y[ns];
						}
					}
					assert(Y.count(s));
					model.addConstr(Y[s] <= expr);
				}
			}
		}
	}
	// cerr << Z.size() << " " << 11 * 11 * 6 * 13 << endl;

	for (int t = maxT-1; t >= 0; t--) {
		cerr << t << endl;
		for (int x = minx; x <= maxx; x++) {
			for (int y = miny; y <= maxy; y++) {
				// // No two robots can collide, unless they are moving in the same direction
				if (t < maxT) {
					for (int d0 = 1; d0 <= 4; d0++) {
						assert(Z.count(node(x, y, t, d0)));
						GRBLinExpr cons = Z[node(x, y, t, d0)];
						for (int d1 = 0; d1 <= 4; d1++) {
							if (d0 == d1) continue;
							if (!in_bounds(x + dx[d0], y + dy[d0])) continue;
							
							assert(Z.count(node(x + dx[d0], y + dy[d0], t, d1)));
							cons += Z[node(x + dx[d0], y + dy[d0], t, d1)];
						}
						model.addConstr(cons <= 1);
					}
				}

				// Each robot must pick a direction
				for (int r = 0; r < ins.n; r++) {
					GRBLinExpr expr;
					for (int d = 0; d <= 4; d++) {
						state u(r, x, y, t, d);
						assert(Y.count(u));
						expr += Y[u];
					}
					state v(r, x, y, t, 5);
					assert(Y.count(v));
					// stringstream ss;
					// ss << "The value at " << v.n.x << " " << v.n.y << " " << v.n.t << " must diffuse";
					model.addConstr(expr == Y[v]);//, ss.str());
				}
			}
		}
	}

	// Each robot must go somewhere
	for (int t = maxT-1; t >= 0; t--) {
		for (int i = 0; i < ins.n; i++) {
			vector<GRBVar> terms;
			vector<double> coeffs;
			for (int x = minx; x <= maxx; x++) {
				for (int y = miny; y <= maxy; y++) {
					terms.push_back(Y[state(i, x, y, t, 5)]);
					coeffs.push_back(1);				
				}
			}

			GRBLinExpr expr;
			expr.addTerms(&coeffs[0], &terms[0], terms.size());
			model.addConstr(expr == 1);//, "Each robot must go somewhere");
		}
	}

	cerr << "Done adding variables" << endl;

	// input output objectives
	for (int i = 0; i < ins.n; i++) {
		pt s = ins.start[i];
		pt t = ins.target[i];
		state st(i, s.x, s.y, 0, 5);
		state en(i, t.x, t.y, maxT, 5);
		// cerr << "Start " << i << ": " << s.x << " " << s.y << endl;
		// cerr << "End " << i << ": " << t.x << " " << t.y << endl;
		assert(Y.count(st) && Y.count(en));
		model.addConstr( Y[st] == 1);//, "start is 1" );
		model.addConstr( Y[en] == 1);//, "start is 1" );
	}

	for (int i = 0; i < ins.n; i++) {
		for (int x = minx; x <= maxx; x++) {
			for (int y = miny; y <= maxy; y++) {
				if (x != ins.start[i].x || y != ins.start[i].y) {
					state st(i, x, y, 0, 5);
					assert(Y.count(st));
					model.addConstr( Y[st] == 0 );//, "non-start is 0" );
				}

				if (x != ins.target[i].x || y != ins.target[i].y) {
					state en(i, x, y, maxT, 5);
					assert(Y.count(en));
					model.addConstr( Y[en] == 0 );//, "non-start is 0");
				}
			}
		}
	}

	// // makespan minimizing objective
	// vector<GRBVar> d(maxT+1);
	// cerr << d.size() << endl;
	// for (int t = maxT; t >= 0; t--) {
	// 	d[t] = model.addVar(0, 1, 1, GRB_BINARY);//, "d" + to_string(t));
	// 	GRBLinExpr expr;
	// 	for (int r = 0; r < ins.n; r++) {
	// 		for (int x = minx; x <= maxx; x++) {
	// 			for (int y = miny; y <= maxy; y++) {
	// 				for (int d = 1; d <= 4; d++) {
	// 					expr += Y[state(r, x, y, t, d)];
	// 				}
	// 			}
	// 		}
	// 	}
	// 	model.addConstr( ins.n * 30 * 30 * 4 * d[t] >= expr );
	// }

	cerr << "About to optimize..." << endl;
	model.set(GRB_IntParam_SolutionLimit, 1);
	model.set(GRB_IntParam_MIPFocus, 1);
	// model.set(GRB_DoubleParam_TimeLimit, 60);
	model.optimize();
	// model.computeIIS();
	// model.write("test.ilp");
	
	// for (auto kv : Y) {
	// 	auto s = kv.first;
	// 	if (kv.second.get(GRB_DoubleAttr_X) > 0) {
	// 		cerr << s.r << " " << s.n.t << " " << s.n.x << " " << s.n.y << " " << s.n.d << endl;
	// 	}
	// }

	// Let's trace the paths
	ins.time = maxT;
	ins.moves.resize(maxT);
	vector<pt> c = ins.start;
	for (int t = 0; t < maxT; t++) {
		ins.moves[t].resize(ins.n);
		for (int i = 0; i < ins.n; i++) {
			//cerr << "looking for move at time " << t << " for robot " << i << endl;
			//cerr << "current position: " << c[i].x << " " << c[i].y << endl;
			int d = 0;
			while (true) {
				assert(Y.count(state(i, c[i].x, c[i].y, t, d)));
				auto& var = Y[state(i, c[i].x, c[i].y, t, d)];
				if (var.get(GRB_DoubleAttr_X) > 0) {
					break;
				}
				d++;
			}
			c[i] = c[i] + dxy[d];
			ins.moves[t][i] = d;
		}
	}
	ins.write();
	return 0;
}
back to top