Hydrothermal Scheduling Problem with Conic Programming

Description: Hydrothermal Scheduling Problem using Second-Order Cones

Tags: amplpy, conic, second-order cone, quadratic cone, nonlinear programming, scheduling, engineering, power generation, geothermal energy, hydropower

Notebook author: Gleb Belov <gleb@ampl.com>

References: Wood, A J, and Wollenberg, B F, Example Problem 7b. In Power Generation, Operation and Control. John Wiley and Sons, 1984, p. 202.

# Install dependencies
%pip install -q amplpy pandas

# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["coin", "mosek"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

Hydrothermal scheduling problem involves allocating the total power demand and losses among the hydro and thermal generators in a least-cost way. The scheduling period is typically a few days long. The hydraulic flow constraints and the limits on generator outputs have to be observed in the scheduling problem.

AMPL model

%%writefile hydrothermal.mod

param nperiods;        ## Number of periods (12 hours long)

set TT := 0..nperiods;
set T  := 1..nperiods;

param load {T};        ## load (MW) for the t-th period
param losscof;         ## loss coeff for hydro generation / 0.00008 /
param nhours           ## number of hours in each period  / 12      /
    default 12;

param vLU {1..2};      ## storage volume range
param thermalLU {1..2};## steam output range
param hydroUB;         ## hydro output upper bound

var cost_sqr{T}>=0;    ## quadratic term of the costs
var thermal{T}         ## output from the steam thermal plant (MW)
               >=thermalLU[1], <=thermalLU[2];
var hydro{T}           ## output from the hydro plant         (MW)
               >=0, <=hydroUB;
var loss{T}    >=0;    ## total loss                          (MW)
var q{TT}      >=0;    ## hydro flow rate in acre-ft per hour
var v{TT}              ## reservoir storage volume at the end of t
               >=vLU[1], <=vLU[2];

minimize Cost:
    1.15*nhours*nperiods * sum {t in T} (500 + 8*thermal[t] + cost_sqr[t]);
    
s.t. CostPerPeriodSqr {t in T}:  ## Extract quadratic terms of the objective
    cost_sqr[t] >= 0.0016*(thermal[t] ** 2); ## into a second-order cone

s.t. Loss {t in T}:    ## loss calculated as function of hydro output
    loss[t] >= losscof*(hydro[t] ** 2);      ## another cone

s.t. Demand {t in T}:  ## demand plus loss must be met from hydro and thermal
    thermal[t] + hydro[t] == load[t] + loss[t];

s.t. Flow {t in T}:    ## hydraulic continuity equation
    v[t] == v[t-1] + (2000 - q[t]) * nhours;

s.t. Dischar {t in T}: ## calculation of hydro discharge
    q[t] == 330 +4.97*hydro[t];

Overwriting hydrothermal.mod

Load data directly from Python data structures using amplpy

m = AMPL()
m.read("hydrothermal.mod")

m.param["nperiods"] = 6

m.param["load"] = [1200, 1500, 1100, 1800, 950, 1300]
m.param["losscof"] = 0.00008

m.param["vLU"] = [60e3, 120e3]
m.param["thermalLU"] = [150, 1500]
m.param["hydroUB"] = 1000

m.eval("fix v[0] := 100e3;")  ## initial storage volume

Solve with Ipopt

m.option["solver"] = "ipopt"
m.get_output("solve;")
m.obj["Cost"].value()

4366944.151328332

Retrieve solution as a pandas dataframe

m.var["thermal"].to_pandas().T

	1	2	3	4	5	6
thermal.val	839.196498	955.803132	801.935007	1079.933203	734.194	864.004863

m.var["hydro"].to_pandas().T

	1	2	3	4	5	6
hydro.val	371.866263	570.207825	305.533028	767.148093	219.666261	452.365932

Solve with Mosek as a Conic Program

m.option["solver"] = "mosek"
# m.option["mosek_options"] = "outlev=1"
m.solve()

MOSEK 10.0.43: MOSEK 10.0.43: optimal; objective 4366944.122
0 simplex iterations
17 barrier iterations

Solve with Mosek as a Quadratic Program (QCP)

m.option["solver"] = "mosek"
m.option["mosek_options"] = "cvt:socp=0"
m.solve()

MOSEK 10.0.43: cvt:socp=0
MOSEK 10.0.43: optimal; objective 4366944.115
0 simplex iterations
17 barrier iterations