#include <iostream>
#include "gurobi_c++.h"
#include <cmath>
#include <cstdlib>
using namespace std;
//參數定義
const int N = 3;
const int V = 1;
double t[N][N];
double x[V][N][N];
const int TMAX=10000;
double rc[V][N];
double lc[N];
double q[V];
double y[V][N];
double c[V];
double r=5;
double ra[V][N];
double Bmin[V];
double battery[V];
double a=1, b=1, f=1, d=1;
double e[N];
double l[N];
double leave[V][N];
double enter[V][N];
double service[N];
void readfile();
int main()
{
try
{
//Gurobi環境設定
GRBEnv test = GRBEnv();
GRBModel model = GRBModel(test);
readfile();
//決策變數定義
GRBVar x[V][N][N];
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
for (int j = 0;j < N;j++)
{
x[k][i][j]=model.addVar(0, 1, 0, GRB_BINARY);
}
}
}
GRBVar leave[V][N];
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
leave[k][i]=model.addVar(0,GRB_INFINITY,0,GRB_CONTINUOUS);
}
}
GRBVar enter[V][N];
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
enter[k][i]=model.addVar(0, GRB_INFINITY, 0, GRB_CONTINUOUS);
}
}
GRBVar ra[V][N];
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
ra[k][i] = model.addVar(0, GRB_INFINITY, 0, GRB_CONTINUOUS);
}
}
GRBVar y[V][N];
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
y[k][i] = model.addVar(0, GRB_INFINITY, 0, GRB_CONTINUOUS);
}
}
//整合變數與參數進模型
model.update();
//模型限制式撰寫
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
for (int j = 0;j < N;j++)
{
if (i == j)
{
model.addConstr(x[k][i][j] == 0);
}
}
for(int i=0;i<N;i++)
{
GRBLinExpr LHS = 0;
for (int j = 0;j < N;j++)
{
LHS += x[k][i][j];
}
model.addConstr(LHS == 1);
}
for (int j = 0;j < N;j++)
{
GRBLinExpr LHS1 = 0;
for (int i = 0;i < N;i++)
{
LHS1 += x[k][j][i];
}
model.addConstr(LHS1 == 1);
}
for (int i = 0;i < N;i++)
{
GRBLinExpr LHS1 = 0;
GRBLinExpr LHS2 = 0;
for (int j = 0;j < N;j++)
{
LHS1 += x[k][i][j];
LHS2 += x[k][j][i];
}
model.addConstr(LHS1 - LHS2 == 0);
}
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
for (int j = 0;j < N;j++)
{
model.addConstr(leave[k][i] + t[i][j] - enter[k][j] <= TMAX*(1 - x[k][i][j]));
}
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
model.addConstr(leave[k][i] >= enter[k][i] + service[i]);
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
model.addConstr(leave[k][i] >= e[i] + service[i]);
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
model.addConstr(leave[k][i] <= l[i]);
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
model.addConstr(leave[k][i] + t[i][0] <= TMAX);
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
for (int j = 0;j < N;j++)
{
model.addConstr(rc[k][j] <= rc[k][i] - (lc[i] * x[k][i][j]) + (q[k] * (1 - x[k][i][j])));
}
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
model.addConstrs(rc[k][i] <= q[k]);
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
for (int j = 0;j < N;j++)
{
model.addConstr(y[k][i] - (c[k] * t[i][j] * x[k][i][j]) + (battery[k] * (1 - x[k][i][j])) >= y[k][j]);
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
for (int j = 0;j < N;j++)
{
model.addConstr(y[k][i] + ra[k][i] + (battery[k] * (1 - x[k][i][j])) >= y[k][j]);
}
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
model.addConstr(battery[k] - y[k][i] >= ra[k][i]);
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
for (int j = 0;j < N;j++)
{
model.addConstr(enter[k][i] + (r*ra[k][i]) + (t[i][j] * x[k][i][j]) >= enter[k][j] - (TMAX*(1 - x[k][i][j])));
}
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
for (int j = 0;j < N;j++)
{
model.addConstr(Bmin[k] <= leave[k][i] - (c[k] * t[i][j] * x[k][i][j]));
}
}
}
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
model.addConstr(y[k][i] <= battery[k]);
}
}
model.set(GRB_IntAttr_ModelSense, 1);
GRBQuadExpr Obj = 0;
for (int k = 0;k < V;k++)
{
for (int i = 0;i < N;i++)
{
for (int j = 0;j < N;j++)
{
Obj += (a*t[i][j] * x[k][i][j])+(b*(e[i]-enter[k][i])*x[k][i][j])+(f*(r*ra[k][i])*x[k][i][j])+x[k][0][j];
}
}
}
model.setObjective(Obj);
//最佳化模型
model.optimize();
//檢查是否達成最佳化
int status = model.get(GRB_IntAttr_Status);
if (status == GRB_OPTIMAL)
{
//輸出目標式函數以及其值
double ObjValue = model.get(GRB_DoubleAttr_ObjVal);
cout << "目標總成本=" << ObjValue << endl;
}
else if (status == GRB_INF_OR_UNBD) {
cout << "Infeasible or unbounded" << endl;
}
else if (status == GRB_INFEASIBLE) {
cout << "Infeasible" << endl;
}
else if (status == GRB_UNBOUNDED) {
cout << "Uunbounded" << endl;
}
else {
cout << "Optimization was stopped with status" << status << endl;
}
double ObjValue = model.get(GRB_DoubleAttr_ObjVal);
cout << "total cost= " << ObjValue << endl;
double LBound = model.get(GRB_DoubleAttr_ObjBound);
cout << "lower bound= " << LBound << endl;
}
catch (GRBException e) {
cout << "Error code = " << e.getErrorCode() << endl;
cout << e.getMessage() << endl;
}
catch (...) {
cout << "Exception during optimization" << endl;
}
system("PAUSE");
}
void readfile()
{
ifstream fin("traveltime.txt");
ifstream fin1("customerdemand.txt");
ifstream fin2("carloading.txt");
ifstream fin3("carbattery.txt");
ifstream fin4("servicetime.txt");
ifstream fin5("carbatterymin.txt");
ifstream fin6("windowsopen.txt");
ifstream fin7("windowsclose.txt");
ifstream fin8("vehiclecon.txt");
for (int i = 0;i < N;i++)
{
for (int j = 0;j < N;j++)
{
fin >> t[i][j];
//cout << t[i][j]<<" ";
}
//cout<<endl;
}
for (int i = 0;i < V;i++)
{
fin1 >> lc[i];
fin2 >> q[i];
fin3 >> battery[i];
fin4 >> service[i];
fin5 >> Bmin[i];
fin6 >> e[i];
fin7 >> l[i];
fin8 >> c[i];
}
}