forcing a battery model to charge OR discharge

[Patriek Brouwer]

Hi All,

I am trying to force my power in and power out to be one of them positive. Unfortunately it is still not working..

Underneath my constraints given;

def p_in_max(model,t):

    return model.p_in[t] <= (model.ess['ESS']/tau)

model.p_in_max = Constraint(T, rule=p_in_max, doc = 'Maximum power in')

def p_out_max(model,t):

    return model.p_out[t] <= (model.ess['ESS']/tau)

model.p_out_max = Constraint(T, rule=p_out_max, doc = 'Maximum power out')

Is it true that it is not possible in pyomo to say if p_in > 0.0 then p_out = 0.0?

def test_rule(model,t):

    if value(model.p_in[t]> 0.0):

        return model.p_out[t] == 0.0

    elif value(model.p_out[t]< 0.0):

        return model.p_in[t] == 0.0

model.p_out_test = Constraint(T, rule=test_rule, doc = 'testje')

But my solver is still ging an error.

ERROR: evaluating object as numeric value: p_in[0]
        (object: <class 'pyomo.core.base.var._GeneralVarData'>)
    No value for uninitialized NumericValue object p_in[0]
ERROR: evaluating object as numeric value: 0.0  <  p_in[0]
        (object: <class 'pyomo.core.expr.expr_pyomo5.InequalityExpression'>)
    No value for uninitialized NumericValue object p_in[0]
ERROR: Rule failed when generating expression for constraint p_out_test with
    index 0: ValueError: No value for uninitialized NumericValue object
    p_in[0]
ERROR: Constructing component 'p_out_test' from data=None failed: ValueError:
    No value for uninitialized NumericValue object p_in[0]

ValueError: No value for uninitialized NumericValue object p_in[0]

Or is there another way of implementing this in my model?

Thanks in advance!

Best, Patriek

[Vincent Reinbold]

Hi Patriek,

I propose you another formulation of the same problem. Let say we have the following parameters and variable (see the doc for explanations).

m.p      = Var(m.time, doc='energy derivative with respect to time',  initialize=0)

m.pd     = Var(m.time, doc='discharging power',  within=NonNegativeReals, initialize=0)
m.pc     = Var(m.time, doc='charging power',     within=NonNegativeReals, initialize=0)
m.u      = Var(m.time, doc='binary variable',    within=Binary,           initialize=0)
m.pcmax  = Param(default=UB,     doc='maximal charging power',     mutable=True, within=PositiveReals)
m.pdmax  = Param(default=UB,     doc='maximal discharging power',  mutable=True, within=PositiveReals)

def _p_balance(b, t):
    return b.p[t] - b.pd[t] + b.pc[t] == 0

def _pdmax(b, t):
    return b.pd[t] - b.u[t] * b.pdmax <= 0

def _pcmax(b, t):
    return b.pc[t] + b.u[t] * b.pcmax <= b.pcmax

m._pdmax     = Constraint(m.time, rule=_pdmax,     doc='Discharging power bound')
m._pcmax     = Constraint(m.time, rule=_pcmax,     doc='Charging power bound')
m._p_balance = Constraint(m.time, rule=_p_balance, doc='Power balance constraint')

The main idea of this formulation is to use one binary variable 'u' (indexed by the time set) that will describe your conditions (if pd > 0, u = 1 and if pc > 0, u = 0).
You will have to define some efficiency (charging and discharging) and also write a proper energy balance equation. Depending on those value, one could get rid of this binary value. But this is compulsory for the general case. 

Good luck !

Vincent

[Patriek Brouwer]

Hi Vincent,

I am trying this with a binary constraint but it is not working. I don’t understand your explanation. Can you explain?

model.P = Var(T, within = Binary)

M = 10000000

def p_batt1(model,t):

    return model.p_out[t] - model.p_in[t] <= M*(1-model.P[t])

model.p_batt1 = Constraint(T, rule=p_batt1, doc = 'Testje1')

def p_batt2(model,t):

    return model.p_in[t] - model.p_out[t] <= M*model.P[t]

model.p_batt2 = Constraint(T, rule=p_batt2, doc = 'Testje2')

def p_batt3(model,t):

    return eta == eta_in*model.P[t] + eta_out*(1-model.P[t])

model.p_batt3 = Constraint(T, rule=p_batt3, doc = 'Testje3')

-- 
You received this message because you are subscribed to the Google Groups "Pyomo Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pyomo-forum/0298079a-a5ff-4862-9798-b013b201c9cc%40googlegroups.com.

[Vincent Reinbold]

Hi Patriek,

I will try to explain better :

The battery power, named p, has to be described using two positive variables, one for charging 'pc' and one for descharging 'pd'. As a result you need, 3 continuous variables and one binary variable to ensure that pc (and respectively pd), should be 0 when pd (and respectivelly pc) is not null. Thus, p is defined in R, pc and pd in R⁺ and u in {0, 1}. and they are all indexed by the time set noted m.time.

This leads to this definition :

m.p      = Var(m.time, doc='energy derivative with respect to time',  initialize=0)

m.pd     = Var(m.time, doc='discharging power',  within=NonNegativeReals, initialize=0)
m.pc     = Var(m.time, doc='charging power',     within=NonNegativeReals, initialize=0)
m.u      = Var(m.time, doc='binary variable',    within=Binary,           initialize=0)

to ensure that pd = 0 when pc is not, (and the other way arround) we can add the following constraints :

m.pcmax  = Param(default=UB,     doc='maximal charging power',     mutable=True, within=PositiveReals)
m.pdmax  = Param(default=UB,     doc='maximal discharging power',  mutable=True, within=PositiveReals)

def _p_balance(b, t):
    return b.p[t] - b.pd[t] + b.pc[t] == 0

def _pdmax(b, t):
    return b.pd[t] - b.u[t] * b.pdmax <= 0

def _pcmax(b, t):
    return b.pc[t] + b.u[t] * b.pcmax <= b.pcmax

m._pdmax     = Constraint(m.time, rule=_pdmax,     doc='Discharging power bound')
m._pcmax     = Constraint(m.time, rule=_pcmax,     doc='Charging power bound')
m._p_balance = Constraint(m.time, rule=_p_balance, doc='Power balance constraint')

Take a minute to think about those inequality, to be sure it is suits your needs. 

I suppose that, you want to go further and write a proper energy equality constraint, taking into account charging and discharging efficiencies. For that purpose, 
you will need to define, the energy in the battery (as a variable), the efficiencies (as a Param), initial and final energies (as Params), such as: 

m.e      = Var(m.time,            doc='energy in battery',          initialize=0)
m.emin   = Param(default=0,       doc='minimum energy (kWh)',       mutable=True, within=NonNegativeReals)
m.emax   = Param(default=10000,   doc='maximal energy (kWh)',       mutable=True)
m.e0     = Param(default=None,    doc='initial state',              mutable=True)
m.ef     = Param(default=None,    doc='final state',                mutable=True)
m.etac   = Param(default=1.0,     doc='charging efficiency',        mutable=True)
m.etad   = Param(default=1.0,     doc='discharging efficiency',     mutable=True)

def _e_balance(m, t):
    if not (m.etac.value and m.etad.value == 1.):
        return m.de[t] == 1 / 3600 * (m.pc[t] * m.etac - m.pd[t] / m.etad)
    else:
        return m.de[t] == 1 / 3600 * (m.pc[t] - m.pd[t])

def _e_initial(m, t):
    if m.e0.value is not None:
        if t == 0:
            return m.e[t] == m.e0
    return Constraint.Skip

def _e_final(m, t):
    if m.ef.value is None:
        return Constraint.Skip
    if t == m.time.last():
        return m.e[t] == m.ef
    else:
        return Constraint.Skip

def _e_min(m, t):
    if m.emin.value is None:
        return Constraint.Skip
    return m.e[t] >= m.emin

def _e_max(m, t):
    if m.emax.value is None:
        return Constraint.Skip
    return m.e[t] <= m.emax

m._e_balance = Constraint(m.time, rule=_energy_balance, doc='Energy balance constraint')
m._e_initial = Constraint(m.time, rule=_e_initial,      doc='Initial energy constraint')
m._e_final   = Constraint(m.time, rule=_e_final,        doc='Final stored energy constraint')
m._e_min     = Constraint(m.time, rule=_e_min,          doc='Minimal energy constraint')
m._e_max     = Constraint(m.time, rule=_e_max,          doc='Maximal energy constraint')

If this is still not clear to you, I suggest you send your code or minimal problem associated with the details of the errors.

Cheers,

Vincent Reinbold

Le lundi 23 septembre 2019 13:57:36 UTC+2, Patriek Brouwer a écrit :

Hi Vincent,

I am trying this with a binary constraint but it is not working. I don’t understand your explanation. Can you explain?

model.P = Var(T, within = Binary)

M = 10000000

def p_batt1(model,t):

    return model.p_out[t] - model.p_in[t] <= M*(1-model.P[t])

model.p_batt1 = Constraint(T, rule=p_batt1, doc = 'Testje1')

def p_batt2(model,t):

    return model.p_in[t] - model.p_out[t] <= M*model.P[t]

model.p_batt2 = Constraint(T, rule=p_batt2, doc = 'Testje2')

def p_batt3(model,t):

    return eta == eta_in*model.P[t] + eta_out*(1-model.P[t])

model.p_batt3 = Constraint(T, rule=p_batt3, doc = 'Testje3')

To unsubscribe from this group and stop receiving emails from it, send an email to pyomo...@googlegroups.com.

[Patriek Brouwer]

Hi Vincent, 

Many thanks for your help. Already stuck for quite some time now.  Would you have a look at my script below? though I was doing what you suggest.

# Specifying the Solver

opt = solvers.SolverFactory('glpk')

# Declare Concrete Model

model = ConcreteModel()

# loop over input parameters /indexeren

# apx_price_dict = {t:apx_price[t] for t in T}

solar_dict = {t:solar[t] for t in T}

wind_dict = {t:wind[t] for t in T}

# input parameters

model.p_solar = Param(T, initialize=solar_dict, doc = 'Solar power (MW)')

model.p_wind = Param(T, initialize=wind_dict, doc = 'Wind power (MW)')

# variables

model.to_grid = Var(T, bounds=(0.0,export_limit), doc = 'Export (MW)')

model.p_in = Var(T, bounds=(0.0,np.Infinity), doc = 'Charging power (MW)')

model.p_out = Var(T, bounds=(0.0,np.Infinity), doc = 'Discharging power (MW)')

model.soc = Var(T, bounds=(soc_0,np.Infinity), doc = 'State of charge (MWh)')

model.ess = Var(K, bounds=(0.0,np.Infinity), doc = 'Storage size (MWh)')

def soc_max(model,t):

    return model.soc[t] <= model.ess['ESS']

model.soc_max = Constraint(T, rule=soc_max, doc = 'Maximum state-of-charge')

def p_in_max(model,t):

    return model.p_in[t] <= (model.ess['ESS']/tau)

model.p_in_max = Constraint(T, rule=p_in_max, doc = 'Maximum power in')

def p_out_max(model,t):

    return model.p_out[t] <= (model.ess['ESS']/tau)

model.p_out_max = Constraint(T, rule=p_out_max, doc = 'Maximum power out')

model.P = Var(T, within = Binary)

M = 10000000

def p_batt1(model,t):

    return model.p_out[t] - model.p_in[t] <= M*(1-model.P[t])

model.p_batt1 = Constraint(T, rule=p_batt1, doc = 'Testje1')

def p_batt2(model,t):

    return model.p_in[t] - model.p_out[t] <= M*model.P[t]

model.p_batt2 = Constraint(T, rule=p_batt2, doc = 'Testje2')

def p_batt3(model,t):

    return eta == eta_in*model.P[t] + eta_out*(1-model.P[t])

model.p_batt3 = Constraint(T, rule=p_batt3, doc = 'Testje3')

def soc_rule(model,t):

    if(t == 0):

        return model.soc[t] == soc_0 + eta*model.p_in[t] - (1/eta)*model.p_out[t]

    else:

        return model.soc[t] == model.soc[t-1] + eta*model.p_in[t] - (1/eta)*model.p_out[t]

model.state_equation = Constraint(T, rule=soc_rule, doc = 'Battery state of charge equation')

    

def balance_rule(model,t):

    return -model.to_grid[t]- model.p_in[t] + model.p_out[t] + model.p_wind[t] + model.p_solar[t] == 0.0

model.balance_equation = Constraint(T, rule=balance_rule, doc='Load balance')

def objective_rule(model):

    output = - C_annualized_batt*model.ess['ESS'] 

    return output

model.objective = Objective(rule = objective_rule, sense=maximize, doc='Objective function')

def pyomo_postprocess(options=None, instance=None, results=None):

    model.objective.display()

opt = SolverFactory('glpk')

results = opt.solve(model)

results.write()

I think there is something wrong in my p_batt1, p_batt2, or p_batt3 constraints.

To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pyomo-forum/505eabdd-13eb-4ad4-a954-a215a627171f%40googlegroups.com.

[Vincent Reinbold]

Hi Patriek,

Why don't you integrate the variable, parameters and constraints I gave you ?

I would replace the constraint 'p_batt1', 'p_batt2', 'p_batt3', 'p_in_max' and 'p_out_max' by the one I gave you (m._pdmax ,m._pcmax). Also soc_rule does not seem homogeneous  to me, soc is usually normalized (between 0 and 100 %).

Vincent

Le lundi 23 septembre 2019 16:10:12 UTC+2, Patriek Brouwer a écrit :

Hi Vincent, 

To view this discussion on the web visit https://groups.google.com/d/msgid/pyomo-forum/505eabdd-13eb-4ad4-a954-a215a627171f%40googlegroups.com.

[luis gustavo cordero bautista]

Hello Vincent,

I began working on Python, and Pyomo to optimize power flow. I had some experience in AMPL but the reality is that I want to learn pyomo because

I had an interview with a professor who asked what softwares I know. I know Matlab and AMPL but not in great detail. This professor strongly suggest me to start learning Python,

Pyomo, OpenDSS. Unfortunately, I wont have the change for visiting foreign universities this year due to pandemic. However, I am determined to learn from internet, forums ,docs,etc.

I came across this conversation, so, I wonder if you can help with examples. I also have a github where I am uploading my codes. I had a problem translating ampl data to pyomo data,

could give me a hand, this is the code that I translated to pyomo ---> https://github.com/luidis2020/OPF-MV-Distribution-Network

[KRISHAN GOPAL SHARMA]

Hi everyone 

I am also facing the same problem
I am unable to understand this piece of code

[Diarmid Roberts]

Krishnan,

I'd suggest that you try the formulation given by Vincent on 23 September 2019.

Cheers,

Diarmid

To view this discussion on the web visit https://groups.google.com/d/msgid/pyomo-forum/54af3401-0588-4bdf-9f66-279a98ecd130n%40googlegroups.com.

--

Diarmid Roberts

PhD Candidate

The University of Sheffield

Recent Publications:

Flow batteries for energy management: Novel algebraic modelling approaches to properly assess their value