Unit Commitment MINLP with Knitro#
Description: Solving a nonlinear Unit Commitment problem with Knitro using MP features for logic and multi-objective optimization. The goal of this notebook is to show a straightforward and clear way of using nonlinear solvers for complex models with logical expressions and also hierarchical multi-objective optimization.
Tags: mp, knitro, mp2nl, nonlinear, quadratic, minlp, unit-commitment, electric-power-industry, energy, multi-objective, gurobi, xpress, mp2nl
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:
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. 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%.
This will be handled by suffixes for the objectives, .objpriority and .objreltol.
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/down limits
subject to Ramp {gen in GENERATORS, t in TIME: ord(t) > 1}:
ramp_down_limit[gen] <= abs(power_generated[gen,t] - power_generated[gen,prev(t)]) <= ramp_up_limit[gen];
# Objective 1: Minimize total operating + startup cost
minimize 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];
# Objective 2: Minimize total emissions
minimize Total_Emissions:
sum {gen in GENERATORS, t in TIME} emission_rate[gen] * power_generated[gen,t];
# Multi-objective support
suffix objpriority;
suffix objabstol;
suffix objreltol;
Writing 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. See the following example for the knitro solver.
# Hierarchical Multi-objective configuration
# Higher priority first
ampl.eval("let Total_Cost.objpriority := 2;")
# Set 5% degradation tolerance for total cost
ampl.eval("let Total_Cost.objreltol := 0.05;")
# Second objective (less priority)
ampl.eval("let Total_Emissions.objpriority := 1;")
SOLVER = "knitro"
# SOLVER='gurobi'
# SOLVER='xpress'
# Time limit
max_seconds = 30
if SOLVER == "knitro":
ampl.solve(
solver="mp2nl",
knitro_options="maxtime=" + str(max_seconds),
verbose=True,
mp2nl_options="solver=knitro obj:multi=2",
)
else:
ampl.solve(
solver=SOLVER,
mp_options="obj:multi=2 outlev=1 timelim=" + str(max_seconds),
verbose=True,
)
is_committed = ampl.get_data("is_committed").to_pandas()
produce = ampl.get_data("power_generated").to_pandas()
startup = ampl.get_data("is_startup").to_pandas()
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: 7880 | 7880
bounds: lower upper range | lower upper range
0 0 3946 | 0 0 3946
free fixed | free fixed
0 3934 | 0 3934
cont. binary integer | cont. binary integer
4438 3442 0 | 4438 3442 0
Number of constraints: 7432 | 7432
eq. ineq. range | eq. ineq. range
linear: 1505 5927 0 | 1505 5927 0
quadratic: 0 0 0 | 0 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 1008 18638 | 1008 18638
quadratic: 504 0 504 | 504 0 504
total: 1008 18638 504 | 1008 18638 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 290236. 17.9151 36.5787 0.080
2 363883. 2.05722 21.8905 0.091
3 334560. 6.59314e-02 13.8155 0.097
4 253744. 4.97603e-05 3.22554 0.104
5 235374. 1.94157e-05 3.21251 0.110
6 205954. 2.70174e-08 0.979615 0.117
7 195467. 3.61624e-11 0.277754 0.125
8 193454. 6.11788e-12 9.12916e-02 0.133
9 193226. 3.99780e-12 7.38386e-02 0.139
10 192816. 1.60538e-13 2.27965e-02 0.146
11 192769. 1.59872e-14 1.00191e-02 0.153
12 192750. 9.88098e-15 7.85043e-04 0.159
13 192746. 2.66454e-15 1.32173e-04 0.165
14 192746. 6.66134e-16 1.83278e-05 0.173
15 192746. 6.66134e-16 1.83278e-05 0.177
16 192746. 6.66134e-16 1.83278e-05 2.727
Root node cutting planes
------------------------
Iter Cuts Best solution Best bound Gap Time
value value (secs)
---- ---- ------------- ---------- --- ------
0 0 192746. 2.882
1 3564 207566. 11.237
2 3698 216905. 11.609
3 3797 225312. 11.994
4 4160 232576. 12.536
5 4201 237303. 13.087
6 4733 248311. 13.653
7 4962 254138. 14.217
8 5117 258076. 14.690
9 5033 258597. 15.325
10 5099 259186. 15.726
11 5118 259229. 16.136
12 5347 259916. 16.728
13 6028 260817. 17.396
14 5792 260937. 17.847
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 260952. 18.383
1 2 260952. 18.681
1 2 260989. 19.394
11 12 261179. 22.477
20 21 261339. 26.337
37 38 261484. 29.051
47 44 261492. 30.173
EXIT: Time limit reached. No integer feasible point found.
Final Statistics for MIP
------------------------
Final objective value =
Final bound value = 2.61491996798331238e+05
Final optimality gap (abs / rel) = inf / inf
# of root cutting plane rounds = 15
# of restarts = 0
# of nodes processed = 47 (51.923s)
# 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 = 92 (51.923s)
Total program time (secs) = 30.28605 (133.743 CPU time)
Time spent in evaluations (secs) = 0.00000
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 38 / 30
Mixed-integer rounding cuts = 4270 / 4779
Gomory cuts = 0 / 0
Flow-cover cuts = 0 / 0
Probing cuts = 27330 / 11800
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 0 / 0 / 0.000s
Rounding heuristic = 16 / 0 / 0.115s
MPEC heuristic = 0 / 0 / 0.000s
Local search heuristic = 6 / 0 / 0.339s
Fix-and-propagate heuristics = 0 / 0 / 0.000s
===========================================================================
Artelys Knitro 15.1.0: maxtime=30
=======================================
Commercial License
Artelys Knitro 15.1.0
=======================================
Knitro changing mip_method from AUTO to 1.
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: 7880 | 7880
bounds: lower upper range | lower upper range
0 0 3946 | 0 0 3946
free fixed | free fixed
0 3934 | 0 3934
cont. binary integer | cont. binary integer
4438 3442 0 | 4438 3442 0
Number of constraints: 7433 | 7433
eq. ineq. range | eq. ineq. range
linear: 1505 5927 0 | 1505 5927 0
quadratic: 0 1 0 | 0 1 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 432 19646 | 432 19646
quadratic: 0 504 504 | 0 504 504
total: 432 19646 504 | 432 19646 504
Knitro using Branch and Bound method with 8 threads.
Initial points
--------------
Read root node initial point:
Objective value: 3.540231e+02
Feasibility error: 1.896679e+02
Feasible: no
Read MIP primal point:
Objective value: 3.540231e+02
Feasibility error: 1.896679e+02
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, 1e+05] | [1e+00, 2e+04]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
1 708.261 174.246 0.597659 0.070
2 1300.32 149.008 0.498104 0.078
3 1582.94 137.592 0.452925 0.087
4 1825.79 138.338 0.411700 0.095
5 2269.69 526.186 0.424402 0.105
6 2509.89 706.731 0.443513 0.114
7 2700.75 842.311 0.477584 0.122
8 2844.83 959.326 0.516635 0.130
9 2885.46 978.970 0.536805 0.139
10 2903.51 998.351 0.558833 0.147
11 2905.52 1024.41 0.583095 0.156
12 2903.69 1024.77 0.586742 0.164
13 2974.03 1003.03 0.620527 0.174
14 3006.85 993.693 0.631773 0.182
15 3051.16 978.413 0.646000 0.191
16 3068.95 973.076 0.651423 0.199
17 3080.08 969.628 0.654822 0.206
18 3093.65 966.335 0.658525 0.213
19 3108.90 961.940 0.662894 0.219
20 3126.08 957.293 0.667514 0.225
30 3251.81 919.428 0.722830 0.292
40 3337.88 891.037 0.770952 0.354
50 3399.42 872.882 0.792710 0.416
60 3550.16 827.253 0.811430 0.477
70 3703.26 781.200 0.819571 0.535
80 3871.74 731.359 0.824109 0.597
90 4060.21 675.583 0.824864 0.664
100 5173.32 355.712 1.21540 0.745
110 6364.69 1.38311 0.958701 0.825
120 6347.22 7.63297 0.980129 0.927
130 6350.46 12.2919 0.980129 0.983
140 6353.66 17.5784 0.980129 1.041
150 6355.01 20.1849 0.980129 1.101
160 6356.61 21.0418 0.980129 1.162
170 6358.59 20.0746 0.980129 1.222
180 6361.95 18.2051 0.980129 1.280
190 6365.17 15.7874 0.989336 1.337
200 6367.04 14.1843 1.00047 1.394
300 6370.34 0.992993 1.11095 1.943
400 6370.29 0.991425 1.03447 3.417
500 6372.60 1.03841 1.11397 4.018
600 6376.06 1.04003 1.06979 4.658
700 6374.85 2.00466 1.01643 5.239
800 6372.51 2.24411 1.03161 5.839
900 6372.82 2.26277 1.18018 6.351
1000 6372.40 8.09501 1.20258 6.839
1100 6370.63 0.992234 1.17436 7.328
1200 6370.46 1.00900 0.980129 7.889
1300 6370.66 1.27561 1.04526 8.468
1400 6370.89 1.29202 1.08312 9.080
1500 6365.10 2.21506 1.18018 9.587
1600 6363.87 1.28544 1.18018 10.229
1700 6367.95 1.77749 1.18018 10.710
1800 6370.05 0.991950 0.971891 11.160
1900 6370.38 0.988622 0.959296 12.118
2000 6370.93 0.985430 0.981499 12.832
2100 6331.70 17.8134 0.957360 13.384
2200 6370.56 0.992158 1.18018 13.917
2300 6370.50 49727.9 1.00250e+09 14.803
2400 6370.54 99830.1 9.56452e+08 15.352
2500 6378.92 122878. 0.800000 15.981
2600 5360.07 64752.9 2245.04 16.780
2700 3044.68 111.656 14.7455 17.598
2800 3162.45 111.359 0.800000 18.449
2900 5539.09 27.7422 0.800000 19.081
3000 6365.61 90675.0 3.42658e+08 19.824
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 -inf 20.105
3 0 inf 20.975
EXIT: Problem appears to be infeasible.
Final Statistics for MIP
------------------------
Final objective value =
Final bound value = +inf
Final optimality gap (abs / rel) = inf / inf
# of root cutting plane rounds = 0
# of restarts = 0
# of nodes processed = 3 (21.703s)
# 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 = 4 (21.703s)
Total program time (secs) = 20.97487 (21.978 CPU time)
Time spent in evaluations (secs) = 0.00000
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 0 / 0
Mixed-integer rounding cuts = 0 / 0
Gomory cuts = 0 / 0
Flow-cover cuts = 0 / 0
Probing cuts = 0 / 0
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 0 / 0 / 0.000s
Rounding heuristic = 0 / 0 / 0.000s
MPEC heuristic = 0 / 0 / 0.000s
Local search heuristic = 4 / 0 / 0.257s
Fix-and-propagate heuristics = 0 / 0 / 0.000s
===========================================================================
MP2NL 0.1: Knitro 15.1.0: Convergence to an infeasible point. Problem may be locally infeasible.
objective 354.0230613; optimality gap 1.8e+308
3 nodes; 4 subproblem solves
------------ WARNINGS ------------
WARNING: "Tolerance violations"
Type MaxAbs [Name] MaxRel [Name]
* aux var bounds 1E+00 1E+00
* variable integrality 5E-01 -
* aux var integrality 5E-01 -
* expr '_abs' 6E+01 1E+00
* algebraic con(s) 7E-01 -
* algebraic con(s) 2E+02 1E+00
* expr '_and' - -
* expr '_condlineq' - -
* expr '_condlinge' - -
* expr '_not' - -
* expr '_or' - -
*: Using the solver's aux variable values.
Documentation: mp.ampl.com/modeling-tools.html#automatic-solution-check.
suffix numiters OUT;
suffix opterror OUT;
suffix relaxbnd OUT;
suffix feaserror OUT;
suffix incumbent OUT;
suffix numfcevals OUT;
Objective = Total_Cost
print("=== Objective Values ===")
total_cost = ampl.obj["Total_Cost"].value()
total_emissions = ampl.obj["Total_Emissions"].value()
print(f"Total cost: {total_cost:.2f}$")
print(f"Total emissions: {total_emissions:.2f} tons CO₂")
=== Objective Values ===
Total cost: 91867.09$
Total emissions: 354.02 tons CO₂
More documentation of the AMPL/MP interface (in particular, MP2NL).
More energy examples.