Using multiple objectives in your model#
Description: We show how to use multiple objectives with Amplpy using a nonlinear Unit Commitment problem. We won’t be using native or emulated features from the solver interface, but emulating manually a lexicographic multiobjective problem.
Tags: warm-start, mp, multi-objective, gurobi, xpress, knitro, mp2nl, electric-power-industry, unit-commitment
Notebook author: Marcos Dominguez Velad <marcos@ampl.com>
# Install dependencies
%pip install -q amplpy pandas numpy
# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook
ampl = ampl_notebook(
modules=["mp2nl", "knitro", "gurobi", "xpress"], # modules to install
license_uuid="default", # license to use
) # instantiate AMPL object and register magics
The problem#
We are solving a version of Unit Commitment with generators with minimum and maximum outputs and ramp limits. There are linear and quadratic costs, so the problem becomes a MINLP, and we aim to minimize CO2 emissions while minimizing production costs, so it is a multi-objective problem.
We will focus on the nonlinear and multi-objective part of the problem for this case. The model is written in AMPL for it to keep easy to read.
Warning: this notebook uses commercial solvers, so when running on Cloud you may need a license for these solvers (Knitro, Gurobi, Xpress), or run locally.
MINLP#
One of the 2 objectives is minimizing total cost, which is a nonlinear expression.
Total cost: $$ \min ; \sum_{g \in \text{GEN}} \sum_{t \in \text{T}} \left( c_g^{\text{lin}} , p_{g,t}
c_g^{\text{quad}} , p_{g,t}^2 \right) ;+; \sum_{g \in \text{GEN}} \sum_{t \in \text{T}} c_g^{\text{start}} , y_{g,t}$$
Where \(p_{g,t}\) is the amount of energy produced by the generator \(g\) in time \(t\), and \(y_{g,t}\) is 1 if the generator started at period \(t\) (startup cost).
There are also some logical constraints to make the model clear. We have some binary variables to model easier if a generator is turned on or has just started:
Producing if and only if commited (binaries \(x_{g,t}\)): $\( p_{g,t} > 0 \;\Longleftrightarrow\; x_{g,t} = 1 \qquad \forall g \in \text{GEN},\; t \in \text{T} \)$
Commited if and only if producing: $\( y_{g,t} \;\Longleftrightarrow\; \big( x_{g,t} = 1 \;\land\; x_{g,t-1} = 0 \big) \qquad \forall g \in \text{GEN} \)$
Amount of energy produced is either 0 or at least over a threshold (minimum input). $\( p_{g,t} = 0 \;\;\text{or}\;\; p_{g,t} \ge \underline{P}_g \qquad \forall g \in \text{GEN},\; t \in \text{T} \)$
Ramp limits:
\[
R_g^{\downarrow}
\;\le\;
\left| p_{g,t} - p_{g,t-1} \right|
\;\le\;
R_g^{\uparrow}
\qquad \forall g \in \text{GEN}
\]
Solvers don’t necessarily handle logical expressions, but the AMPL/MP interface takes care of reformulating efficiently these expressions.
Multi-objective#
We have 2 objectives, the first is to minimize cost, and the second is to minimize carbon emissions.
We are going to do hierarchical optimization minimizing first the total cost minimization problem, and later the emissions, allowing a degradation of the total cost of 5%.
We will emulate the multi-objective process by adding a constraint that ensure bounds for the already found objectives. The idea is to iterate over the different objectives (2 in this case), solve for a single objective, and add a constraint that ensures that the current objective won’t become worse in the following runs. Then, solve for the next objective.
This can be handled natively with suffixes in the objectives assigning priorities or degradation coefficients. Look at this example for such a version: https://colab.ampl.com/notebooks/unit-commitment-minlp-with-knitro.html
MINLP Unit Commitment model#
%%writefile unit_commitment.mod
set GENERATORS;
set TIME ordered;
param demand {TIME} >= 0; # Power demand at each time
param min_output {GENERATORS} >= 0; # Minimum power output
param max_output {g in GENERATORS} >= min_output[g]; # Maximum power output
param ramp_up_limit {GENERATORS} >= 0;
param ramp_down_limit {GENERATORS} >= 0;
param linear_cost {GENERATORS}; # Linear cost coefficient
param quadratic_cost {GENERATORS} >= 0; # Quadratic cost coefficient
param startup_cost {GENERATORS} >= 0; # Startup cost
param emission_rate {GENERATORS} >= 0; # Tons CO2 per MW produced
var is_committed {GENERATORS, TIME} binary; # 1 if generator is ON
var power_generated {gen in GENERATORS, TIME} >= 0 <= max_output[gen]; # MW produced
var is_startup {GENERATORS, TIME} binary; # 1 if generator starts up
# Generation only if committed
subject to Generation_Commitment {gen in GENERATORS, t in TIME}:
power_generated[gen,t] > 0 <==> is_committed[gen,t] = 1;
# Meet demand in each period
subject to Demand_Satisfaction {t in TIME}:
sum {gen in GENERATORS} power_generated[gen,t] >= demand[t];
# Startup in first period
subject to Startup_First {gen in GENERATORS}:
is_startup[gen, first(TIME)] == is_committed[gen, first(TIME)];
# Startup logic in subsequent periods
subject to Startup_Transition {gen in GENERATORS, t in TIME: ord(t) > 1}:
is_startup[gen,t] <==> (is_committed[gen,t] and !is_committed[gen,prev(t)]);
subject to Min_Gen_If_On {gen in GENERATORS, t in TIME}:
power_generated[gen,t] == 0 or power_generated[gen,t] >= min_output[gen];
# Ramp-up limits
subject to Ramp_Up {gen in GENERATORS, t in TIME: ord(t) > 1}:
power_generated[gen,t] - power_generated[gen,prev(t)] <= ramp_up_limit[gen];
# Ramp-down limits
subject to Ramp_Down {gen in GENERATORS, t in TIME: ord(t) > 1}:
power_generated[gen,prev(t)] - power_generated[gen,t] <= ramp_down_limit[gen];
# Objective 1: Minimize total operating + startup cost
var Total_Cost = sum {gen in GENERATORS, t in TIME}
(linear_cost[gen] * power_generated[gen,t] +
quadratic_cost[gen] * power_generated[gen,t]^2)
+ sum {gen in GENERATORS, t in TIME} startup_cost[gen] * is_startup[gen,t];
minimize Total_Cost_Obj:
Total_Cost;
param Total_Cost_Bound default 1e21;
s.t. Ensure_Total_Cost_Obj:
Total_Cost <= Total_Cost_Bound;
# Objective 2: Minimize total emissions
var Total_Emissions = sum {gen in GENERATORS, t in TIME} emission_rate[gen] * power_generated[gen,t];
minimize Total_Emissions_Obj:
Total_Emissions;
param Total_Emissions_Bound default 1e21;
s.t. Ensure_Total_Emissions_Obj:
Total_Emissions <= Total_Emissions_Bound;
Overwriting unit_commitment.mod
Solving the model#
Generate the problem data#
The AMPL model is isolated from the input data for more readability and maintainance
import pandas as pd
import numpy as np
# Unit Commitment data
generators = ["G1", "G2", "G3", "G4", "G5", "G6", "G7"]
# Generators data
generators_data = pd.DataFrame(
{
"min_output": [20, 30, 25, 15, 10, 40, 0],
"max_output": [100, 120, 90, 60, 50, 150, 30],
"ramp_up_limit": [40, 50, 30, 25, 20, 60, 10],
"ramp_down_limit": [40, 50, 30, 25, 20, 60, 10],
"linear_cost": [20, 16, 18, 22, 24, 14, 12],
"quadratic_cost": [0.04, 0.05, 0.06, 0.03, 0.04, 0.036, 0.1],
"startup_cost": [400, 300, 360, 200, 160, 600, 160],
"emission_rate": [0.7, 0.5, 0.6, 0.4, 0.3, 0.8, 0.0],
},
index=generators,
)
# Generate random demand
num_time_periods = 24 * 3
time_periods = list(range(1, num_time_periods + 1))
np.random.seed(42)
base_demand = 150 + 40 * np.sin(np.linspace(0, 3 * np.pi, num_time_periods))
noise = np.random.normal(0, 10, num_time_periods)
demand = (base_demand + noise).clip(min=100).round().astype(int)
Assign the UC data#
# Optimization model and data
ampl = AMPL()
ampl.read("unit_commitment.mod")
ampl.set["TIME"] = time_periods
ampl.set["GENERATORS"] = generators
ampl.param["demand"] = demand
ampl.set_data(generators_data, "GENERATORS")
Running the solver#
For nonlinear solvers like Knitro, we need the mp2nl interface to make it handle logical constraints and multi-objective.
def run_solver(model, solver, max_seconds, objective):
if solver == "knitro":
model.solve(
solver="mp2nl",
problem=objective,
knitro_options="maxtime=" + str(max_seconds),
verbose=True,
mp2nl_options="solver=knitro obj:multi=2",
)
else:
model.solve(
solver=solver,
problem=objective,
mp_options="obj:multi=2 outlev=1 timelim=" + str(max_seconds),
verbose=True,
)
return model
# We don't need these constraints for the first solve
ampl.con["Ensure_Total_Cost_Obj"].drop()
ampl.con["Ensure_Total_Emissions_Obj"].drop()
ampl = run_solver(
model=ampl, solver="knitro", max_seconds=30, objective="Total_Cost_Obj"
)
# Second objective
# We need to ensure Total Cost keeps low
objective_value = ampl.obj["Total_Cost_Obj"].value()
ampl.param["Total_Cost_Bound"] = objective_value # add the bound from the previous run
ampl.con[
"Ensure_Total_Cost_Obj"
].restore() # restore the constraint that ensures the bound
ampl = run_solver(
model=ampl, solver="knitro", max_seconds=30, objective="Total_Emissions_Obj"
) # minimize Total Emissions
MP2NL 0.1: nl:solver = knitro
obj:multi = 2
Artelys Knitro 15.1.0: maxtime=30
=======================================
Commercial License
Artelys Knitro 15.1.0
=======================================
Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.
concurrent_evals 0
datacheck 0
feastol 1e-06
feastol_abs 1e-06
findiff_numthreads 1
hessian_no_f 1
hessopt 1
maxtime 30
opttol 1e-06
opttol_abs 0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 1.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 0.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_heuristic_fixpropagate from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
Problem Characteristics | Presolved
-----------------------
Problem type: convex MIQP
Objective: minimize / quadratic
Number of variables: 6389 | 6389
bounds: lower upper range | lower upper range
0 0 3449 | 0 0 3449
free fixed | free fixed
0 2940 | 0 2940
cont. binary integer | cont. binary integer
3444 2945 0 | 3444 2945 0
Number of constraints: 5941 | 5941
eq. ineq. range | eq. ineq. range
linear: 1008 4933 0 | 1008 4933 0
quadratic: 0 0 0 | 0 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 1008 13171 | 1008 13171
quadratic: 504 0 504 | 504 0 504
total: 1008 13171 504 | 1008 13171 504
Knitro using Branch and Bound method with 8 threads.
Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.
Coefficient range:
linear objective: [1e+01, 6e+02] | [2e+00, 1e+02]
linear constraints: [1e-04, 2e+02] | [1e-04, 1e+02]
quadratic objective: [3e-02, 1e-01] | [5e-03, 2e-02]
quadratic constraints: [0e+00, 0e+00] | [0e+00, 0e+00]
variable bounds: [1e+00, 2e+02] | [1e+00, 2e+02]
constraint bounds: [1e+00, 2e+02] | [1e+00, 2e+02]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
1 291358. 8.57743 37.0858 0.030
2 310343. 1.28140 14.5466 0.036
3 294269. 3.18702e-02 17.1532 0.040
4 249437. 6.63367e-05 8.72878 0.046
5 223518. 2.96053e-06 3.43428 0.051
6 201867. 6.69924e-09 0.872545 0.057
7 194258. 1.58015e-11 0.195604 0.063
8 193366. 4.52656e-12 0.136435 0.069
9 193067. 2.14607e-12 0.103684 0.074
10 192834. 3.75433e-13 5.90346e-02 0.079
11 192772. 2.56328e-14 8.37065e-03 0.084
12 192753. 1.26565e-14 2.55204e-03 0.090
13 192747. 3.55271e-15 2.36541e-04 0.095
14 192746. 6.66134e-16 2.50623e-05 0.100
15 192746. 6.66134e-16 2.50623e-05 0.103
16 192746. 6.66134e-16 2.50623e-05 0.596
Root node cutting planes
------------------------
Iter Cuts Best solution Best bound Gap Time
value value (secs)
---- ---- ------------- ---------- --- ------
0 0 258322. LS 192746. 25.39% 0.668
1 1373 258322. 193014. 25.28% 0.844
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 258322. 193019. 25.28% 10.553
1 2 243020. LS 193019. 20.57% 10.777
6 7 243020. 193064. 20.56% 11.578
23 24 243020. 193215. 20.49% 13.197
46 47 243020. 193274. 20.47% 14.952
68 69 243020. 193280. 20.47% 16.337
94 95 243020. 193280. 20.47% 18.010
124 125 243020. 193280. 20.47% 19.448
146 147 243020. 193280. 20.47% 21.087
175 176 243020. 193280. 20.47% 22.619
204 205 243020. 193280. 20.47% 24.411
237 238 243020. 193285. 20.47% 26.198
267 268 243020. 193286. 20.46% 27.563
308 309 243020. 193295. 20.46% 29.026
339 340 243020. 193388. 20.42% 30.002
EXIT: Time limit reached. Integer feasible point found.
Final Statistics for MIP
------------------------
Final objective value = 2.43019557358517312e+05
Final bound value = 1.93388030559908279e+05
Final optimality gap (abs / rel) = 49631.5 / 0.204229 (20.42%)
# of root cutting plane rounds = 2
# of restarts = 0
# of nodes processed = 339 (90.742s)
# of strong branching evaluations = 0 (0.000s)
# of function evaluations = 0 (0.000s)
# of gradient evaluations = 0 (0.000s)
# of hessian evaluations = 0 (0.000s)
# of hessian-vector evaluations = 0
# of subproblems processed = 787 (127.562s)
Total program time (secs) = 30.02850 (136.621 CPU time)
Time spent in evaluations (secs) = 0.00000
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 72 / 28
Mixed-integer rounding cuts = 7 / 7
Gomory cuts = 0 / 0
Flow-cover cuts = 0 / 0
Probing cuts = 4811 / 1711
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 2 / 0 / 29.724s
Rounding heuristic = 3 / 0 / 0.139s
MPEC heuristic = 2 / 0 / 7.878s
Local search heuristic = 6 / 5 / 0.085s
Fix-and-propagate heuristics = 0 / 0 / 0.000s
===========================================================================
MP2NL 0.1: Knitro 15.1.0: Time limit reached. Current point is feasible.
objective 243019.5574; optimality gap 4.96e+04
339 nodes; 587 subproblem solves
suffix numiters OUT;
suffix opterror OUT;
suffix relaxbnd OUT;
suffix feaserror OUT;
suffix incumbent OUT;
suffix numfcevals OUT;
MP2NL 0.1: nl:solver = knitro
obj:multi = 2
Artelys Knitro 15.1.0: maxtime=30
=======================================
Commercial License
Artelys Knitro 15.1.0
=======================================
Knitro changing mip_method from AUTO to 1.
MINLP solver shifted start point to satisfy bounds (1434 variables).
concurrent_evals 0
datacheck 0
feastol 1e-06
feastol_abs 1e-06
findiff_numthreads 1
hessian_no_f 1
hessopt 1
maxtime 30
opttol 1e-06
opttol_abs 0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 1.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 0.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_heuristic_fixpropagate from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
Problem Characteristics | Presolved
-----------------------
Problem type: convex MIQCQP
Objective: minimize / linear
Number of variables: 6389 | 6389
bounds: lower upper range | lower upper range
0 0 3449 | 0 0 3449
free fixed | free fixed
0 2940 | 0 2940
cont. binary integer | cont. binary integer
3444 2945 0 | 3444 2945 0
Number of constraints: 5942 | 5942
eq. ineq. range | eq. ineq. range
linear: 1008 4933 0 | 1008 4933 0
quadratic: 0 1 0 | 0 1 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 432 14179 | 432 14179
quadratic: 0 504 504 | 0 504 504
total: 432 14179 504 | 432 14179 504
Knitro using Branch and Bound method with 8 threads.
Initial points
--------------
Read root node initial point:
Objective value: 7.500004e+03
Feasibility error: 1.000000e+00
Feasible: no
Read MIP primal point:
Objective value: 7.500004e+03
Feasibility error: 1.000000e+00
Feasible: no, Knitro will try to repair it
Coefficient range:
linear objective: [3e-01, 8e-01] | [3e-01, 8e-01]
linear constraints: [1e-04, 6e+02] | [1e-04, 1e+02]
quadratic objective: [0e+00, 0e+00] | [0e+00, 0e+00]
quadratic constraints: [3e-02, 1e-01] | [5e-03, 2e-02]
variable bounds: [1e+00, 2e+02] | [1e+00, 2e+02]
constraint bounds: [1e+00, 2e+05] | [1e+00, 4e+04]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
1 7448.40 1.00269 0.800000 0.088
2 7084.28 382.169 0.800000 0.098
3 6868.54 268.304 0.800000 0.109
4 6223.27 221.301 0.800000 0.120
5 5941.09 344.059 0.779046 0.130
6 5697.69 219.374 0.655226 0.141
7 5234.44 191.858 0.468373 0.150
8 4524.04 5.03593e-05 0.317117 0.158
9 3825.39 2.22045e-16 0.171111 0.166
10 3607.84 2.22045e-16 9.44585e-02 0.179
11 3590.21 2.22045e-16 9.09358e-02 0.187
12 3579.93 2.22045e-16 4.90865e-03 0.197
13 3575.03 2.22045e-16 6.10587e-03 0.205
14 3572.97 2.22045e-16 6.28583e-04 0.214
15 3572.61 2.22045e-16 9.31141e-04 0.223
16 3572.14 2.22045e-16 2.00080e-04 0.233
17 3571.97 9.54547e-04 1.91392e-05 0.241
18 3571.95 2.22045e-16 4.35289e-06 0.250
19 3571.95 2.22045e-16 3.75991e-07 0.257
20 3571.95 2.22045e-16 3.75991e-07 0.262
Root node cutting planes
------------------------
Iter Cuts Best solution Best bound Gap Time
value value (secs)
---- ---- ------------- ---------- --- ------
0 0 7335.80 LS 3571.95 51.31% 0.347
1 1352 7335.80 3577.75 51.23% 0.668
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 7335.80 3577.75 51.23% 0.900
1 2 5460.50 RFP 3577.75 34.48% 1.313
3 4 5460.50 3577.75 34.48% 2.571
17 18 5460.50 3577.75 34.48% 6.204
31 32 5460.50 3579.04 34.46% 9.266
43 42 5460.50 3579.94 34.44% 11.959
70 65 5460.50 3582.11 34.40% 15.191
100 95 5460.50 3582.39 34.39% 18.555
138 133 5460.50 3582.57 34.39% 21.288
170 163 5460.50 3583.33 34.38% 23.957
177 170 3585.49 FCRD 3584.00 0.04% 25.093
188 179 3584.78 FCRD 3584.00 0.02% 26.321
197 176 3584.78 3584.00 0.02% 26.590
222 159 3584.78 3584.43 0.01% 27.889
EXIT: Optimal solution found.
Final Statistics for MIP
------------------------
Final objective value = 3.58478331934105518e+03
Final bound value = 3.58442851101700535e+03
Final optimality gap (abs / rel) = 0.354808 / 9.89762e-05 (0.01%)
# of root cutting plane rounds = 2
# of restarts = 0
# of nodes processed = 222 (103.505s)
# of strong branching evaluations = 0 (0.000s)
# of function evaluations = 0 (0.000s)
# of gradient evaluations = 0 (0.000s)
# of hessian evaluations = 0 (0.000s)
# of hessian-vector evaluations = 0
# of subproblems processed = 262 (126.237s)
Total program time (secs) = 28.32510 (184.969 CPU time)
Time spent in evaluations (secs) = 0.00000
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 40 / 16
Mixed-integer rounding cuts = 33 / 48
Gomory cuts = 0 / 0
Flow-cover cuts = 0 / 0
Probing cuts = 4580 / 5618
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 1 / 1 / 0.165s
Rounding heuristic = 5 / 2 / 0.153s
MPEC heuristic = 2 / 0 / 22.530s
Local search heuristic = 6 / 1 / 0.395s
Fix-and-propagate heuristics = 0 / 0 / 0.000s
===========================================================================
MP2NL 0.1: Knitro 15.1.0: Locally optimal or satisfactory solution.
objective 3584.783319; optimality gap 0.355
222 nodes; 262 subproblem solves
print("=== Objective Values ===")
total_cost = ampl.obj["Total_Cost_Obj"].value()
total_emissions = ampl.obj["Total_Emissions_Obj"].value()
print(f"Total cost: {total_cost:.2f}$")
print(f"Total emissions: {total_emissions:.2f} tons CO₂")
=== Objective Values ===
Total cost: 243019.56$
Total emissions: 3584.78 tons CO₂
More documentation of the AMPL/MP interface (in particular, MP2NL).
More energy examples.