Customzie Formulation¶

Encapsuled the optimization problem calss, AMS provides direct access to the optimization formulation, where users have the option to customize the formulation without playing with the source code.

In this example, we will walk through the optimization problem structure and show how to customize the formulation.

In [1]:

import numpy as np

import ams

import datetime

In [2]:

print("Last run time:", datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S"))

print(f'ams:{ams.__version__}')

Last run time: 2024-06-18 20:17:45
ams:0.9.8

In [3]:

ams.config_logger(stream_level=20)

Inspect the Optimization Problem Structure¶

In [4]:

sp = ams.load(ams.get_case('5bus/pjm5bus_demo.xlsx'),
              setup=True,
              no_output=True,)

Parsing input file "/Users/jinningwang/Documents/work/mambaforge/envs/amsre/lib/python3.9/site-packages/ams/cases/5bus/pjm5bus_demo.xlsx"...
Input file parsed in 0.0920 seconds.
Zero line rates detacted in rate_b, rate_c, adjusted to 999.
System set up in 0.0021 seconds.

In AMS, a routine collects the descriptive dispatch formulations. DCOPF, RTED, etc, are the subclasses of RoutineBase.

In [5]:

sp.DCOPF.init()

<DCOPF> initialized in 0.0080 seconds.

Out[5]:

True

After successful initialization, the attribute om is populated with CVXPY-based optimization problem.

The user can even hack to the source prob attribute to customize it if necessary.

In [6]:

type(sp.DCOPF.om.prob)

Out[6]:

cvxpy.problems.problem.Problem

Customize Built-in Formulation¶

Here we extend DCOPF with consideration of CO2 emission, where the original formulation can be found in the documentation Routine Reference - DCOPF. To simplify the demonstration, following assumptions are made:

Variable $\boxed{e_g}$ is the CO2 emission of each generator. It is proportional to the generation, described by a parameter $\boxed{k_e}$ in the unit t/p.u..
Total CO2 emmission is limited by a constant cap $\boxed{t_e}$, in the unit $t$.
A tax $\boxed{c_e}$ is imposed on each unit of CO2 emission in the unit of $/p.u., and the tax is included in the objective function.

Thus, the revised formulation is as follows, where box indicates the revision:

min. $\sum ( c_{2} p_{g}^2 + c_{1} p_{g} + u_{g} c_{0} + \boxed{c_{e} e_{g}} )$

s.t.

$\boxed{ e_{g} - k_{e} p_{g} = 0}$

$\boxed{ \sum e_{g} - t_{e} \leq 0}$

$-p_g + c_{trl,ne}p_{g,0} + c_{trl,e}p_{g,\min} \leq 0$

$p_g - c_{trl,ne}p_{g,0} - c_{trl,e}p_{g,\max} \leq 0$

$B_{bus}\theta_{bus} + p^{inj}_{bus} + C_{lpd} + C_{sh}g_{sh} - C_{p}g_{p} = 0$

$-B_f\theta_{bus} - p^{inj}_f - R_{ATEA} \leq 0$

$B_f\theta_{bus} + p^{inj}_f - R_{ATEA} \leq 0$

$-C^T_f\theta_{bus} - \theta_{\max} \leq 0$

$C^T_f\theta_{bus} - \theta_{\max} \leq 0$

Decision variables: $p_g$, $\theta_{bus}$, $\boxed{e_g}$

Note that line flow $p_{lf}$ is calculated as $B_f\theta_{bus} + p^{inj}_f$ after solving the problem.

Add services¶

Services are used to store values or build matrix for easier formulation.

In [7]:

sp.DCOPF.addService(name='te', tex_name='t_e',
                    unit='t', info='emission cap',
                    value=12)

Out[7]:

ValueService: DCOPF.te

Add parameters¶

We need the following parameters to be defined as RParam: ke and ce. They should be 1D array in the same length as the number of generators and te is a scalar.

For a general RParam, it has attributes model, indexer, and imodel to describe its source model and index model. The definition of c2 in DCOPF source code is a good example. However, for ones defined through API, since there is no model containing it, all above attributes are not applicable, and the user should be aware of the sequence of the parameters.

Considering the sequence can be indexed by the generator index, it is used to reference the variables order. Assuming ke is reciprocal to the generator capacity, and ce is the same for each generator, we can define the parameters as follows:

In [8]:

# get the generator indices
stg_idx = sp.DCOPF.pg.get_idx()

# get the value of pmax
pmax = sp.DCOPF.get(src='pmax', attr='v', idx=stg_idx)

# assume the emission factor is 1 for all generators
ke = np.ones_like(pmax)

# assume tax is reciprocal of pmax
ce = np.reciprocal(pmax)

In [9]:

sp.DCOPF.addRParam(name='ke', tex_name='k_e',
                   info='gen emission factor',
                   model=None, src=None, unit=None,
                   v=ke)

Out[9]:

RParam: DCOPF.ke

In [10]:

sp.DCOPF.addRParam(name='ce', tex_name='c_e',
                   info='gen emission tax',
                   model=None, src=None, unit=None,
                   v=ce)

Out[10]:

RParam: DCOPF.ce

Add variables¶

The gerator emission eg is added as a new variable.

In [11]:

sp.DCOPF.addVars(name='eg', tex_name='e_g',
                 info='Gen emission', unit='t',
                 model='StaticGen', src=None)

Out[11]:

Var: StaticGen.eg

Add constraints¶

The CO2 emission is an equality constraint, and the CO2 emission cap is a simple linear inequality constraint.

If wish to revise an existing built-in constraint, you can redefine the constraint e_str attribute.

In [12]:

sp.DCOPF.addConstrs(name='egb', info='Gen emission balance',
                    e_str='eg - mul(ke, pg)', is_eq=True)

Out[12]:

Constraint: egb [ON]

In [13]:

sp.DCOPF.addConstrs(name='eub', info='emission upper bound',
                    e_str='sum(eg) - te', is_eq=False,)

Out[13]:

Constraint: eub [ON]

Revise the objective function¶

The e_str can be revised to include the CO2 emission tax. Here we only need to append the tax term to the original objective function.

In [14]:

sp.DCOPF.obj.e_str += '+ sum(mul(ce, pg))'

Finalize the Customization¶

After revising the problem, remember to initialize it before solving.

In [15]:

sp.DCOPF.init()

<DCOPF> initialized in 0.0101 seconds.

Out[15]:

True

Solve it and Check the Results¶

In [16]:

sp.DCOPF.run(solver='CLARABEL')

<DCOPF> solved as optimal in 0.0102 seconds, converged in 8 iterations with CLARABEL.

Out[16]:

True

Inspect the results.

In [17]:

sp.DCOPF.eg.v

Out[17]:

array([0.20000047, 0.92389198, 5.2       , 3.67610756])

In [18]:

sp.DCOPF.pg.v

Out[18]:

array([0.20000047, 0.92389198, 5.2       , 3.67610756])

In [19]:

sp.DCOPF.obj.v

Out[19]:

5.337040792868603

Load the original problem as a baseline for comparison.

In [20]:

sp0 = ams.load(ams.get_case('5bus/pjm5bus_demo.xlsx'),
               setup=True,
               no_output=True,)

Parsing input file "/Users/jinningwang/Documents/work/mambaforge/envs/amsre/lib/python3.9/site-packages/ams/cases/5bus/pjm5bus_demo.xlsx"...
Input file parsed in 0.0342 seconds.
Zero line rates detacted in rate_b, rate_c, adjusted to 999.
System set up in 0.0021 seconds.

In [21]:

sp0.DCOPF.run(solver='CLARABEL')

<DCOPF> initialized in 0.0058 seconds.
<DCOPF> solved as optimal in 0.0083 seconds, converged in 8 iterations with CLARABEL.

Out[21]:

True

From the comparasion, we can see that the generation schedule changes.

In [22]:

sp0.DCOPF.pg.v

Out[22]:

array([2. , 2.1, 5.2, 0.7])

In [23]:

sp0.DCOPF.obj.v

Out[23]:

2.3445000001347758