The C++ implementation uses IBM CPLEX for optimization. Below is a detailed breakdown of the main code components.
#include <ilcplex/ilocplex.h>
#include <chrono>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
ILOSTLBEGIN
typedef IloArray<IloNumVarArray> NumVar2D;
struct DieselGenerator {
double cost_per_mwh; // Cost ($/MWh)
double p_min; // Minimum capacity (MW)
double p_max; // Maximum capacity (MW)
};
int main(int, char**)
{
// The main function content is below
}
const int T = 24; // Time horizon (hours)
const int numBuses = 4, numDGs = 2, numLoads = 2, numBESS = 1, numRES = 1;
const double price_scale = 0.8;
const bool consider_opf = false;
const double x = 0.1, r = 0.0;
const double line_capacity = 500.0;
std::vector<std::pair<int, int>> lines = {
{0, 1}, {1, 2}, {2, 3}, {3, 0}, {0, 2}
};
// Susceptance and B matrix
std::vector<std::vector<double>> susceptance(numBuses, std::vector<double>(numBuses, 0.0));
std::vector<std::vector<double>> B(numBuses, std::vector<double>(numBuses, 0.0));
// Fill susceptance matrix
for (const auto& line : lines) {
int i = line.first, j = line.second;
double b_ij = -x / (r * r + x * x);
susceptance[i][j] = b_ij;
susceptance[j][i] = b_ij;
}
// Fill B matrix (DC power flow)
for (int i = 0; i < numBuses; ++i) {
for (int j = 0; j < numBuses; ++j) {
if (i != j) B[i][j] = -susceptance[i][j];
else {
double sum = 0.0;
for (int k = 0; k < numBuses; ++k)
if (k != i) sum += susceptance[i][k];
B[i][j] = sum;
}
}
}// Initialize location vectors
std::vector<int> grid_location(numBuses, 0);
std::vector<int> battery_location(numBuses, 0);
std::vector<int> res_location(numBuses, 0);
std::vector<std::vector<int>> dg_locations(numDGs, std::vector<int>(numBuses, 0));
std::vector<std::vector<int>> load_locations(numLoads, std::vector<int>(numBuses, 0));
// Set locations
grid_location[3] = 1;
battery_location[0] = 1;
res_location[0] = 1;
dg_locations[0][1] = 1;
dg_locations[1][2] = 1;
load_locations[0][1] = 1;
load_locations[1][2] = 1;
// Electricity purchase price from the grid ($/MWh) for each hour
std::vector<double> c_buy = {
90, 90, 90, 90, 90, 90,
110, 110, 110, 110, 110, 125,
125, 125, 125, 125, 125, 125,
80, 80, 80, 80, 80, 80
};
// Load demand (MW) per load per hour
std::vector<std::vector<double>> p_load(numLoads, std::vector<double>{
169, 175, 179, 171, 181, 172,
270, 264, 273, 281, 193, 158,
161, 162, 250, 260, 267, 271,
284, 167, 128, 134, 144, 150
});
// Renewable energy forecast (MW) per hour
std::vector<double> p_res = {
16, 17, 17, 100, 181, 172,
270, 264, 273, 281, 193, 158,
161, 162, 250, 260, 267, 271,
284, 167, 128, 134, 144, 150
};
// Battery Energy Storage System (BESS)
const double p_bess_max = 200.0; // Max charge/discharge power (MW)
const double eta = 0.95; // Round-trip efficiency
const double soc_min = 0.10; // Minimum state of charge (fraction)
const double soc_max = 0.90; // Maximum state of charge (fraction)
// Diesel Generator Parameters
std::vector<DieselGenerator> dgs = {
{80.0, 50.0, 200.0}, // DG1: cost, min, max
{70.0, 50.0, 200.0} // DG2: cost, min, max
};
NumVar2D p_dg(env, numDGs);
for (int i = 0; i < numDGs; ++i) {
p_dg[i] = IloNumVarArray(env, T, 0, IloInfinity, ILOFLOAT);
}
IloNumVarArray p_buy(env, T, 0, IloInfinity, ILOFLOAT);
IloNumVarArray p_sell(env, T, 0, IloInfinity, ILOFLOAT);
IloNumVarArray p_bess_chg(env, T, 0, IloInfinity, ILOFLOAT);
IloNumVarArray p_bess_dch(env, T, 0, IloInfinity, ILOFLOAT);
IloNumVarArray soc(env, T, 0, 1, ILOFLOAT);
IloNumVarArray omega(env, T, 0, IloInfinity, ILOINT);
NumVar2D theta(env, numBuses);
for (int i = 0; i < numBuses; ++i) {
theta[i] = IloNumVarArray(env, T, -IloInfinity, IloInfinity, ILOFLOAT);
}IloExpr objective(env);
for (int t = 0; t < T; t++) {
objective += c_buy[t] * p_buy[t] - price_scale * c_buy[t] * p_sell[t];
for (int i = 0; i < numDGs; ++i) {
objective += dgs[i].cost_per_mwh * p_dg[i][t];
}
}
model.add(IloMinimize(env, objective));
for (int t = 0; t < T; t++) {
if (consider_opf) {
// Energy balance with DC power flow
for (int i = 0; i < numBuses; i++) {
IloExpr flow(env);
for (int j = 0; j < numBuses; j++) {
flow += B[i][j] * (theta[i][t] - theta[j][t]);
}
model.add(
(p_buy[t] - p_sell[t]) * grid_location[i]
+ (p_bess_dch[t] - p_bess_chg[t]) * battery_location[i]
+ p_res[t] * res_location[i]
+ p_dg[0][t] * dg_locations[0][i]
+ p_dg[1][t] * dg_locations[1][i]
- p_load[0][t] * load_locations[0][i]
- p_load[1][t] * load_locations[1][i]
- flow == 0
);
}
// Line limits
for (int i = 0; i < numBuses; ++i) {
for (int j = 0; j < numBuses; ++j) {
model.add(B[i][j] * (theta[i][t] - theta[j][t]) <= line_capacity);
model.add(B[i][j] * (theta[i][t] - theta[j][t]) >= -line_capacity);
}
}
// Theta limits
for (int i = 0; i < numBuses; ++i) {
model.add(theta[i][t] >= theta_min);
model.add(theta[i][t] <= theta_max);
}
} else {
// Energy balance constraint without DC power flow
model.add(
p_buy[t] + p_res[t] + p_dg[0][t] + p_dg[1][t] + p_bess_dch[t]
== p_bess_chg[t] + p_load[0][t] + p_load[1][t] + p_sell[t]
);
}
// BESS constraints
if (t == 0) {
model.add(soc[t] == soc[T - 1] + (p_bess_chg[t] * eta - p_bess_dch[t] / eta) / p_bess_max);
model.add(p_bess_chg[t] <= p_bess_max * (1 - soc[T - 1]) / eta);
model.add(p_bess_dch[t] <= p_bess_max * soc[T - 1] * eta);
} else {
model.add(soc[t] == soc[t - 1] + (p_bess_chg[t] * eta - p_bess_dch[t] / eta) / p_bess_max);
model.add(p_bess_chg[t] <= p_bess_max * (1 - soc[t - 1]) / eta);
model.add(p_bess_dch[t] <= p_bess_max * soc[t - 1] * eta);
}
// Binary variable constraints for mutually exclusive charge/discharge
model.add(p_bess_chg[t] <= omega[t] * M);
model.add(p_bess_dch[t] <= (1 - omega[t]) * M);
// State of Charge (SoC) limits
model.add(soc[t] >= soc_min);
model.add(soc[t] <= soc_max);
// Diesel Generator constraints
for (int i = 0; i < numDGs; ++i) {
model.add(p_dg[i][t] >= dgs[i].p_min);
model.add(p_dg[i][t] <= dgs[i].p_max);
}
}
IloCplex cplex(env);
cplex.extract(model);
// cplex.setOut(env.getNullStream());
if (!cplex.solve()) {
env.error() << "Failed" << endl;
throw(-1);
}
double obj = cplex.getObjValue();
std::cout << "Minimized Objective Function: " << obj << std::endl;