CODEWITHSUNDEEP
pythonmathematical-optimizationpyomo
errenmike1806
errenmike1806

Reputation: 109

Redefine variable domain in pyomo

I am a little bit confused about how to use instance to change the model in pyomo. For example, I have a model then I create a variable

    #non negative variable
    SS.lambda_DA_E=pe.Var(SS.T, domain=pe.PositiveReals)

then I create an instance and solve the model

    instance=SS.create_instance()
    solver = po.SolverFactory('mpec_nlp')
    results = solver.solve(instance, load_solutions=True, tee=True)

Now I want to redefine the variable as a free variable, assign it the initial values that I got from the previous solution and solve it again, I tried to use

    def _init_rule_L_E2(m, t):
        return (pe.value(instance.lambda_DA_E[t]))
    instance.lambda_DA_E=pe.Var(SS.T, domain=pe.Reals,  initialize=_init_rule_L_E2)

and then I got an error ValueError: No value for uninitialized NumericValue object lambda_DA_E[t1]

Do I redefine variable correctly? Because I also got a warning

WARNING: Implicitly replacing the Component attribute lambda_DA_E (type=<class
    'pyomo.core.base.var.IndexedVar'>) on block unknown with a new Component
    (type=<class 'pyomo.core.base.var.IndexedVar'>). This is usually
    indicative of a modelling error. To avoid this warning, use
    block.del_component() and block.add_component().

Here is the full code ```

    import pandas as pd
    import pyomo.environ as pe
    import pyomo.opt as po
    import pyomo.mpec as mpc
    
    #DATA
    T=2;
    G=2;
    
    time = ['t{0}'.format(t+1) for t in range(T)]
    gen=['G{0}'.format(g+1) for g in range(G)]
    
    
    #Technical characteristic
    
    elec_maxprod= {'G1': 200, 'G2': 200,} #Only for Generators
    
    D_E={'t1':350, 't2': 150} #El demand
    
    C_fuel = {'G1': 15.0, 'G2': 50.0}
    #MODEL
    SS=pe.ConcreteModel()
    #Duals are desired
    SS.dual = pe.Suffix(direction=pe.Suffix.IMPORT) 
    
    
    ### SETS
    SS.T = pe.Set(initialize = time)
    SS.G = pe.Set(initialize = gen ) 
    
    ### PARAMETERS
    
    SS.P_max = pe.Param(SS.G, initialize = elec_maxprod) #Only for Generators
    SS.D_E = pe.Param(SS.T, initialize = D_E) #El demand
    SS.C_fuel = pe.Param(SS.G, initialize = C_fuel)
    
    ### VARIABLES
    
    SS.p_DA=pe.Var(SS.G, SS.T, domain=pe.PositiveReals)
    
    ##Dual variables
    #positive variables
    
    SS.mu_min_G=pe.Var(SS.G, SS.T, domain=pe.PositiveReals)
    SS.mu_max_G=pe.Var(SS.G, SS.T, domain=pe.PositiveReals)
    
    #free variables
    SS.lambda_DA_E=pe.Var(SS.T, domain=pe.PositiveReals) # day-ahead electricity price in period t [$ per MWh]
    
    ##CONSTRAINTS
    
    ##Maximum and Minimum Generation 
    SS.comp1 = mpc.ComplementarityList()
    for t in SS.T:
        for g in SS.G:
           SS.comp1.add(expr=mpc.complements( SS.p_DA[g,t] >= 0, SS.mu_min_G[g,t] >= 0))
            
    SS.comp2 = mpc.ComplementarityList()
    for t in SS.T:
        for g in SS.G:
            SS.comp2.add(expr=mpc.complements( SS.P_max[g] - SS.p_DA[g,t] >= 0, SS.mu_max_G[g,t] >= 0))
            
    ##Stationarity         
    SS.L1_p_DA = pe.ConstraintList()
    for t in SS.T:
        for g in SS.G:
            SS.L1_p_DA.add( SS.C_fuel[g] - SS.mu_min_G[g,t] + SS.mu_max_G[g,t] - SS.lambda_DA_E[t]== 0)       
            
    #El Balance
    SS.El_DA_bal=pe.ConstraintList()
    for t in SS.T:
        SS.El_DA_bal.add( sum(SS.p_DA[g,t] for g in SS.G) - SS.D_E[t]  == 0 )
            
    
    instance=SS.create_instance() 
    solver = po.SolverFactory('mpec_nlp')
    results = solver.solve(instance, load_solutions=True, tee=True)
    
    print("\nDisplaying Solution\n" + '-'*60)     
    ##Printing the variables results
    for t in SS.T:
        print('price_DA[',t,']', pe.value(instance.lambda_DA_E[t]))
    
    ##Redefine variable with a new feasible region and new initial value    
    def _init_rule_L_E2(m, t):
        return (pe.value(instance.lambda_DA_E[t]))
    instance.lambda_DA_E=pe.Var(SS.T, domain=pe.Reals,  initialize=_init_rule_L_E2) # day-ahead electricity price in period t [$ per MWh]
    
    ##Solve it again
    results = solver.solve(instance, load_solutions=True, tee=True)
    ##Printing the variables results
    for t in SS.T:
        print('price_DA[',t,']', pe.value(instance.lambda_DA_E[t]))

Upvotes: 0

Views: 810

Answers (1)

pybegginer
pybegginer

Reputation: 525

The error is because you're not modifying the current component lambda_DA_E but repleacing it with a new one.

If you want to change any attribute of the component, you may use the assignement like this:

#After solving the whole model.
>>>instance.lambda_DA_E[t].display()
lambda_DA_E : Size=2, Index=T
Key : Lower : Value : Upper : Fixed : Stale : Domain
 t1 :     0 :  None :  None : False :  True : PositiveReals
 t2 :     0 :  None :  None : False :  True : PositiveReals
#After solving the first time
>>>print('changing Domain:')
>>>instance.lambda_DA_E[t].domain=pe.Reals
>>>results = solver.solve(instance, load_solutions=True, tee=True)
>>>instance.lambda_DA_E[t].display()
lambda_DA_E : Size=2, Index=T
Key : Lower : Value              : Upper : Fixed : Stale : Domain
 t1 :  None :  50.00000000097778 :  None : False : False :  Reals
 t2 :  None : 15.000000000977778 :  None : False : False :  Reals

Then you can solve the model again with solver.solve(instance, ...)

Just a couple of annotations:

  1. Since SS is a ConcreteModel() you cannot call instance=SS.create_instance(). Depending on the pyomo version, this may raise a WARNING:

     WARNING: DEPRECATED: Cannot call Model.create_instance() on a constructed model; returning a clone of the current model instance.

    If you want to create a copy of the SS ConcreteModel you can use clone():

     instance = SS.clone()

  2. I have never worked with ComplementarityList from pyomo.environ.mpc, but trying to solve this problem will return me a Warning:

     WARNING: Implicitly replacing the Component attribute v (type=<class 'pyomo.core.base.var.SimpleVar'>) on block comp2[2] with a new Component (type=<class 'pyomo.core.base.var.SimpleVar'>). This is usually indicative of a modelling error. To avoid this warning, use block.del_component() and block.add_component().

I think that mpec_nlp solver use some sort of reformulation that invoked again when resolving. This may not be a problem, but FYI.

Upvotes: 2

Related Questions