Implementation of optimization routine

Support
1

Hi everyone, relative newbie here. Trying to implement the attached algo with the following code:

IMG_1585
IMG_15851920×1330 419 KB

from amplpy import AMPL, ampl_notebook

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

%%ampl_eval

%%ampl_eval

define parameters

param R_cap;      # Rated power capacity of the peaker in MWh
param D;            # Duration of the BESS in hours
param eta;  # Round-trip efficiency of the BESS
param n;
param P_batt;            # BESS price per kW
param z{i in 1..n};      # State of charge at time i in kWh
param Q{i in 1..n};      # Discharging power at time i in kW
param P_elec{i in 1..n}; # Electricity price at time i in $/kWh

%%ampl_eval

define variables

var R;                    # (2) Power capacity of the BESS in kW
var r{i in 1..n} ;        # Charging power at time i in kW

%%ampl_eval

%%ampl_eval

define model and constraints

minimize Total_Cost:
P_batt*R+sum{i in 1..n}P_elec[i]*r[i]; # Minimize total cost of BESS and electricity (1)

subject to charging_constraint{i in 1..n}: r[i]>=0;
subject to charging_constraint2{i in 1..n}: r[i]<=R;       #(2)
subject to max_cap_constraint: R<=R_cap;            # (2)
subject to stat_of_charge{i in 1..n}: z[i]>=0;
subject to stat_of_charge2{i in 1..n}: z[i]<=R*D;           # (4)
subject to charge_discharge_balance{i in 2..n}:
z[i]=z[i-1]+eta*r[i-1]-(1/eta)*Q[i-1];             # (3)
subject to charge_discharge{i in 1..n}: r[i]*Q[i] = 0;        # (6)
subject to discharge_constraint{i in 1..n}: R>=Q[i];        # (7)

# load parameter data
ampl.param["R_cap"]=498
ampl.param["P_batt"]=400
ampl.param["D"]=4
ampl.param["n"]=24 
ampl.param["eta"]=(0.85)**(1/2)
ampl.param["z[1]"]=[0.5*498*4,0] 
ampl.param["P_elec"]=[32.01,41.6,39.73,41.53,39.44,36.6,32.64,48.51,40.86,47.52,62.18,58.29,55.09,60.33,58.34,63.97,67.81,67.5,54.49,50.7,56.87,42.58,49.49,51.18]
ampl.param["Q"]=[8.14,8.14,8.14,10.72,10.16,8.14,8.14,35.86,20.18,15.33,14.79,13.84,12.35,28.77,48.03,46.81,56.72,124.83,140.30,140.16,139.89,128.36,97.42,60.52]

Getting the following error:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[97], line 7
      5 ampl.param["n"]=24 
      6 ampl.param["eta"]=(0.85)**(1/2)
----> 7 ampl.param["z[1]"]=[0.5*498*4,0] 
      8 ampl.param["P_elec"]=[32.01,41.6,39.73,41.53,39.44,36.6,32.64,48.51,40.86,47.52,62.18,58.29,55.09,60.33,58.34,63.97,67.81,67.5,54.49,50.7,56.87,42.58,49.49,51.18]
      9 ampl.param["Q"]=[8.14,8.14,8.14,10.72,10.16,8.14,8.14,35.86,20.18,15.33,14.79,13.84,12.35,28.77,48.03,46.81,56.72,124.83,140.30,140.16,139.89,128.36,97.42,60.52] 

File amplpy/ampl.pyx:883, in amplpy.ampl.AMPL._param.Parameters.__setitem__()

File amplpy/ampl.pyx:320, in amplpy.ampl.AMPL.get_parameter()

File amplpy/util.pxi:55, in amplpy.ampl.PY_AMPL_CALL()

KeyError: 'An entity called z[1] cannot be found.'
2

Here are the conditions that platform wouldn’t allow me to post on first newbie post:

IMG_1586
IMG_15861920×2222 700 KB

3

To assign 0.5*498*4 to z[1], you can write:

ampl.param["z"][1]=0.5*498*4

Or if you also want to assign 0 to z[2] through z[24]. then you can use:

ampl.param["z"]=[0.5*498*4]+[0]*23
4

You were blocked by an automatic spam detector that is built into Discourse. The block has been removed now, and you should be able to reply.

5

Bob, thanks. That helps.

Interestingly, still getting the all zeros solution for R and r. But that violates the constraint:

subject to discharge_constraint{i in 1..n}: R>=Q[i]; # (7)

Thanks, Gui

6

Hi Gui,

Before the solver is called, variables are all zero (unless you have assigned them some other values). So, first check that the solver is being called — by ampl.solve(solver="highs") — before you view the values of R and r in the solution.

If that is not the problem, check the result message from the solver. It should report an optimal solution, as in this example:

HiGHS 1.11.0: optimal solution; objective 88.2

If the result message shows anything else, the solve may have failed, leaving the variables at 0. If you need help, post the complete result message.

If HiGHS does report an optimal solution, then there may be an error in the model or in the data. You can post the complete amplpy program to get some more suggestions.

7

Hi Bob,

Indeed, getting an ‘infeasible’ solution…

Trying to see if I can get a data set that works…

Thanks,
G

8

Hi Bob,

This is a simpler routine:

I’m trying to maximize the profit of a battery that can charge r_i or discharge q_i, with wholesale electricity prices P_elec,i, within a n=24 hour period. (eta is just the efficiency of the battery, and z, the state of charge, which should cap q).

I’ve setup a very simple, do-it-by-hand test case with very high and very low P_elec, so it should be obvious how to maximize profit. But it still only returns 0s.

Any insights highly appreciated!

Thanks, Gui

image
image1191×719 38.3 KB

9

image
image1186×681 31.5 KB

10

The result of display z; shows that parameter z[i] is zero for all i in 1..n. Thus the constraint 0<=q[i]<=z[i] forces all of the variables q[i] to zero; and from that it follows that the optimal values of all the variables r[i] are zero as well.

You’ll need to study the model and your amplpy program to determine how to fix this. As a start, you could consider whether the data for z is specified incorrectly, or whether z ought to be a variable instead of a parameter.

Also, the HiGHS solver does not recognize quadratic constraints, so it is not going to be able to handle r[i]*q[i] = 0 except by making a piecewise-linear approximation. To get an exact result, you can write the constraint instead as

subject to charge_discharge {i in 1..n}:
   r[i] = 0 or q[i] = 0;

AMPL’s solver interface converts this to an equivalent linear constraint, so that HiGHS can handle it.