Assigmnent of negative/positive value of a variable to another variable not working

[Cord Kaldemeyer]

In my simplified model I have a (continous) variable with a lower bound LB below zero and an upper bound UB above zero. Now I want to assign the variable value to other variables depending on the value that the variable has taken.

In Pseudocode the logic can be expressed as follows:

LB > 0
UB > 0
-LB <= Variable1 <= UB

if Variable1 => 0:
    Variable2 = Variable1
    Variable3 = 0
else:
    Variable2 = 0
    Variable3 = abs(Variable1)

Using unequations it could be expressed like this with the use of a binary variable:

LB > 0
UB > 0

-LB <= Variable1 <= UB

# Binary variable that equals 1 if Variable1 > 0 and 0 if Variable1 < 0
Variable1 <= UB * BinaryVar
LB * (1 - BinaryVar) <= Variable1

# Value assignment
Variable2 = Variable1 * BinaryVar
Variable3 = -Variable1 * (1 - BinaryVar)

This of course leads to a quadratic term which can be linearized later..

But if I implement this logic in my model, the value assigment somehow does not work although the binary variable gets set rightly in all cases.

Here's a minimal example with random numbers that should work with glpk or gurobi:

Code # -*- coding: utf-8 -*-

# IMPORTS #####################################################################


from pyomo.environ import *
from pyomo.opt import SolverFactory
import numpy as np
import random
import time


# MODEL CREATION & SETS #######################################################


# Create model
mdl = ConcreteModel()

# Timesteps
timesteps_max = 24
timesteps = range(0, timesteps_max)

# Set for timesteps
mdl.T = Set(ordered=True, initialize=timesteps)


# PARAMETERS ##################################################################


## Technical ------------------------------------------------------------------

# Turbine: Installed capacity in MW
p_el_turb_max = 321

# Compressor: Installed capacity in MW
p_el_comp_max = 68


## Market data (random numbers) -----------------------------------------------


# Positive tertiary reserve: Working prices in EUR/MWh
c_price_el_tr_pos_wrk = dict(zip(timesteps, [random.randrange(0, 450, 1)
                                 for _ in range(timesteps_max)]))
mdl.c_price_el_tr_pos_wrk = Param(mdl.T, initialize=c_price_el_tr_pos_wrk)


# Positive tertiary reserve: Activation level in MW
p_el_tr_pos_act = dict(zip(timesteps, [random.randrange(0, p_el_turb_max, 1)
                           for _ in range(timesteps_max)]))
mdl.p_el_tr_pos_act = Param(mdl.T, initialize=p_el_tr_pos_act)


# Negative tertiary reserve: Working prices in EUR/MWh
c_price_el_tr_neg_wrk = dict(zip(timesteps, [random.randrange(0, 250, 1)
                                 for _ in range(timesteps_max)]))
mdl.c_price_el_tr_neg_wrk = Param(mdl.T, initialize=c_price_el_tr_neg_wrk)


# Spot market: Price in EUR/MWh
c_price_el_spot = dict(zip(timesteps,
                           [random.randrange(-250, 250, 1) for _ in
                            range(timesteps_max)]))
mdl.c_price_el_spot = Param(mdl.T, initialize=c_price_el_spot)


# VARIABLES ###################################################################


# Physical power plant capacity in MW
mdl.p_el_pp_phy = Var(mdl.T, within=Reals, bounds=(-p_el_comp_max,
                                                   p_el_turb_max))

# Capacity at spot markets in MW
mdl.p_el_spot = Var(mdl.T, within=Reals, bounds=(-p_el_comp_max,
                                                 p_el_turb_max))

# Turbine capacity in MW
mdl.p_el_turb = Var(mdl.T, within=Reals)

# Compressor capacity in MW
mdl.p_el_comp = Var(mdl.T, within=Reals)

# Binary variable to assign the turbine and compressor capacity
mdl.p_el_turb_comp_bin = Var(mdl.T, within=Binary)

# Positive tertiary reserve: Provided capacity in MW
mdl.p_el_tr_pos_wrk = Var(mdl.T, within=NonNegativeReals,
                          bounds=(0, p_el_turb_max + p_el_comp_max))

# Negative tertiary reserve: Provided capacity in MW
mdl.p_el_tr_neg_wrk = Var(mdl.T, within=NonNegativeReals,
                          bounds=(0, p_el_comp_max + p_el_turb_max))


# CONSTRAINTS #################################################################


# Physical power plant capacity
def p_el_pp_phy_constr_rule(mdl, t):
    return (mdl.p_el_pp_phy[t] == mdl.p_el_spot[t] + mdl.p_el_tr_pos_wrk[t] -
            mdl.p_el_tr_neg_wrk[t])
mdl.p_el_pp_phy_constr = Constraint(mdl.T, rule=p_el_pp_phy_constr_rule)


# Turbine and compressor capacity definition
def p_el_turb_comp_constr_0_rule(mdl, t):
    return (mdl.p_el_pp_phy[t] <= p_el_turb_max * mdl.p_el_turb_comp_bin[t])
mdl.p_el_turb_comp_constr_0 = Constraint(mdl.T,
                                         rule=p_el_turb_comp_constr_0_rule)


def p_el_turb_comp_constr_1_rule(mdl, t):
    return (-p_el_comp_max * (1 - mdl.p_el_turb_comp_bin[t]) <=
            mdl.p_el_pp_phy[t])
mdl.p_el_turb_comp_constr_1 = Constraint(mdl.T,
                                         rule=p_el_turb_comp_constr_1_rule)


def p_el_turb_comp_constr_2_rule(mdl, t):
    return (mdl.p_el_turb[t] == mdl.p_el_pp_phy[t] * mdl.p_el_turb_comp_bin[t])
mdl.p_el_turb_comp_constr_2 = Constraint(mdl.T,
                                         rule=p_el_turb_comp_constr_2_rule)


def p_el_turb_comp_constr_3_rule(mdl, t):
    return (mdl.p_el_comp[t] == -mdl.p_el_pp_phy[t] *
            (1 - mdl.p_el_turb_comp_bin[t]))
mdl.p_el_turb_comp_constr_3 = Constraint(mdl.T,
                                         rule=p_el_turb_comp_constr_3_rule)


# OBJECTIVES ##################################################################

def obj_rule(mdl):
    return sum(mdl.p_el_spot[t] * mdl.c_price_el_spot[t] -
               mdl.p_el_tr_neg_wrk[t] * mdl.c_price_el_tr_neg_wrk[t] -
               mdl.p_el_tr_pos_wrk[t] * mdl.c_price_el_tr_pos_wrk[t]
               for t in mdl.T)
mdl.obj = Objective(sense=minimize, rule=obj_rule)


# INSTANCE CREATION ###########################################################


# Take the time
tmp_start_time = time.time()

# Create instance and print model

instance = mdl.create()

print("Model creation: %g seconds" % (time.time() - tmp_start_time))

# Choose solver (standard = "glpk")
opt = SolverFactory("gurobi")

# Set options
opt.options["threads"] = 3  # Works only with gurobi
opt.options["mipgap"] = 0.01

results = opt.solve(instance, tee=True, warmstart=True)

print("Solver: %g seconds" % (time.time() - tmp_start_time))


# POSTPROCESSING ##############################################################


# Emtpy dictionary to store results
result_dict = dict()

# Load results in model instance and store them in a dictionary
instance.load(results)
for v in instance.active_components(Var):
    varobject = getattr(instance, v)
    tmp_lst = []
    for index in varobject:
        tmp_lst.append(varobject[index].value)
    result_dict[v] = [tmp_lst]

# Print results
print result_dict

In my Model Variable1 is named "mdl.p_el_pp_phy[t]" and the binary var "mdl.p_el_turb_comp_bin[t]" gets set rightly. The variables "mdl.p_el_turb[t]" and "mdl.p_el_comp[t]" equal Variable2 and Variable3 from the example. Somehow they are always set to zero instead of the value via the binary var...


And I have tried so many different things but cannot figure out why...

Can anyone help? I got really stuck here..

Thanks in advance!
Cord

[David Woodruff]

I don't understand why a model with quadratic constraints should work with gurobi or glpk. You text has the words
"which can be linearized later" but I don't see that you have done that linearization.

On Mon, Feb 16, 2015 at 10:17 AM, Cord Kaldemeyer <cordm...@gmail.com> wrote:

In my simplified model I have a (continous) variable with a lower bound LB below zero and an upper bound UB above zero. Now I want to assign the variable value to other variables depending on the value that the variable has taken.

In Pseudocode the logic can be expressed as follows:

LB > 0
UB > 0
-LB <= Variable1 <= UB

if Variable1 => 0:
    Variable2 = Variable1
    Variable3 = 0
else:
    Variable2 = 0
    Variable3 = abs(Variable1)

Code hier eingeben...

--
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.
For more options, visit https://groups.google.com/d/optout.

[Cord Kaldemeyer]

Sorry, my mistake. The problem cannot be solved by glpk but gurobi can solve quadratic problems. I could also post the linearized version of the problem. Do you think this would make a difference?

Cord

...

[Cord Kaldemeyer]

I have created a linearized version and now it is working somehow! Thanks anyway for your answer which brought me indirectly to my solution ;-)

Cheers
Cord

Linearized version of the model:

# -*- coding: utf-8 -*-

# IMPORTS #####################################################################


from pyomo.environ import *
from pyomo.opt import SolverFactory
import numpy as np
import random
import time


# MODEL CREATION & SETS #######################################################


# Create model
mdl = ConcreteModel()

# Timesteps

timesteps_max = 24 * 7

timesteps = range(0, timesteps_max)

# Set for timesteps
mdl.T = Set(ordered=True, initialize=timesteps)


# PARAMETERS ##################################################################


## Technical ------------------------------------------------------------------

# Turbine: Installed capacity in MW
p_el_turb_max = 321

# Compressor: Installed capacity in MW
p_el_comp_max = 68


## Market data (random numbers) -----------------------------------------------


# Positive tertiary reserve: Working prices in EUR/MWh

c_price_el_tr_pos_wrk = dict(zip(timesteps, [random.randrange(0, 100, 1)

                                 for _ in range(timesteps_max)]))
mdl.c_price_el_tr_pos_wrk = Param(mdl.T, initialize=c_price_el_tr_pos_wrk)


# Positive tertiary reserve: Activation level in MW
p_el_tr_pos_act = dict(zip(timesteps, [random.randrange(0, p_el_turb_max, 1)
                           for _ in range(timesteps_max)]))
mdl.p_el_tr_pos_act = Param(mdl.T, initialize=p_el_tr_pos_act)


# Negative tertiary reserve: Working prices in EUR/MWh

c_price_el_tr_neg_wrk = dict(zip(timesteps, [random.randrange(0, 100, 1)

# Linearization of the quadratic term
# mdl.p_el_turb[t] == mdl.p_el_pp_phy[t] * mdl.p_el_turb_comp_bin[t]
# see http://orinanobworld.blogspot.de/2010/10/
# binary-variables-and-quadratic-terms.html)
def p_el_turb_comp_constr_2_rule(mdl, t):
    return (mdl.p_el_turb[t] <= p_el_turb_max * mdl.p_el_turb_comp_bin[t])

mdl.p_el_turb_comp_constr_2 = Constraint(mdl.T,
                                         rule=p_el_turb_comp_constr_2_rule)


def p_el_turb_comp_constr_3_rule(mdl, t):

    return (mdl.p_el_turb[t] >= -p_el_comp_max * mdl.p_el_turb_comp_bin[t])
mdl.p_el_turb_comp_constr_3 = Constraint(mdl.T,
                                         rule=p_el_turb_comp_constr_3_rule)


def p_el_turb_comp_constr_4_rule(mdl, t):
    return (mdl.p_el_turb[t] <= mdl.p_el_pp_phy[t] - (-p_el_comp_max) *

            (1 - mdl.p_el_turb_comp_bin[t]))

mdl.p_el_turb_comp_constr_4 = Constraint(mdl.T,
                                         rule=p_el_turb_comp_constr_4_rule)


def p_el_turb_comp_constr_5_rule(mdl, t):
    return (mdl.p_el_turb[t] >= mdl.p_el_pp_phy[t] - p_el_turb_max *

            (1 - mdl.p_el_turb_comp_bin[t]))

mdl.p_el_turb_comp_constr_5 = Constraint(mdl.T,
                                         rule=p_el_turb_comp_constr_5_rule)


def p_el_turb_comp_constr_6_rule(mdl, t):
    return (mdl.p_el_comp[t] == mdl.p_el_turb[t] - mdl.p_el_pp_phy[t])
mdl.p_el_turb_comp_constr_6 = Constraint(mdl.T,
                                         rule=p_el_turb_comp_constr_6_rule)

# OBJECTIVES ##################################################################

def obj_rule(mdl):
    return sum(mdl.p_el_spot[t] * mdl.c_price_el_spot[t] -
               mdl.p_el_tr_neg_wrk[t] * mdl.c_price_el_tr_neg_wrk[t] -
               mdl.p_el_tr_pos_wrk[t] * mdl.c_price_el_tr_pos_wrk[t]
               for t in mdl.T)
mdl.obj = Objective(sense=minimize, rule=obj_rule)


# INSTANCE CREATION ###########################################################


# Take the time
tmp_start_time = time.time()

# Create instance and print model
instance = mdl.create()

print("Model creation: %g seconds" % (time.time() - tmp_start_time))

# Choose solver (standard = "glpk")

opt = SolverFactory("glpk")

results = opt.solve(instance)

print("Solver: %g seconds" % (time.time() - tmp_start_time))


# POSTPROCESSING ##############################################################


# Emtpy dictionary to store results

result_dict = dict()

# Load results in model instance and store them in a dictionary
instance.load(results)
for v in instance.active_components(Var):
    varobject = getattr(instance, v)
    tmp_lst = []
    for index in varobject:
        tmp_lst.append(varobject[index].value)
    result_dict[v] = [tmp_lst]

# Print results
print result_dict

...