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.

import numpy as np

import ams

import datetime

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

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

Last run time: 2024-04-21 17:30:31
ams:0.9.6.post3.dev0+gfd3786e

ams.config_logger(stream_level=20)

Inspect the Optimization Problem Structure

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

Parsing input file "/Users/jinningwang/Documents/work/ams/ams/cases/5bus/pjm5bus_demo.xlsx"...
Input file parsed in 0.1127 seconds.
Zero line rates detacted in rate_b, rate_c, adjusted to 999.
System set up in 0.0020 seconds.

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

sp.DCOPF.init()

<DCOPF> initialized in 0.0093 seconds.

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.

type(sp.DCOPF.om.prob)

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: 1. 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.. 1. Total CO2 emmission is limited by a constant cap \(\boxed{t_e}\), in the unit \(t\). 1. 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.

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

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:

# 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)

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

RParam: DCOPF.ke

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

RParam: DCOPF.ce

Add variables

The gerator emission eg is added as a new variable.

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

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.

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

Constraint: egb [ON]

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

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.

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

Finalize the Customization

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

sp.DCOPF.init()

<DCOPF> initialized in 0.0071 seconds.

True

Solve it and Check the Results

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

<DCOPF> solved as optimal in 0.0116 seconds, converged in 10 iterations with ECOS.

True

Inspect the results.

sp.DCOPF.eg.v

array([0.20000037, 0.92389203, 5.2       , 3.6761076 ])

sp.DCOPF.pg.v

array([0.20000037, 0.92389203, 5.2       , 3.6761076 ])

sp.DCOPF.obj.v

5.337040792067892

Load the original problem as a baseline for comparison.

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

Parsing input file "/Users/jinningwang/Documents/work/ams/ams/cases/5bus/pjm5bus_demo.xlsx"...
Input file parsed in 0.0385 seconds.
Zero line rates detacted in rate_b, rate_c, adjusted to 999.
System set up in 0.0023 seconds.

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

<DCOPF> initialized in 0.0071 seconds.
<DCOPF> solved as optimal in 0.0091 seconds, converged in 10 iterations with ECOS.

True

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

sp0.DCOPF.pg.v

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

sp0.DCOPF.obj.v

2.3445000000435057