Python pyomo.environ.Constraint() Examples

def operationMode2(self, pyM, esM, constrName, constrSetName, opVarName, opRateName='operationRateFix',
                       isStateOfCharge=False):
        """
        Define operation mode 2. The operation [commodityUnit*h] is equal to the installed capacity multiplied
        with a time series in:\n
        * [commodityUnit*h] (for storages) or in
        * [commodityUnit] multiplied by the hours per time step (else).\n
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        opVar, capVar = getattr(pyM, opVarName + '_' + abbrvName), getattr(pyM, 'cap_' + abbrvName)
        constrSet2 = getattr(pyM, constrSetName + '2_' + abbrvName)
        factor = 1 if isStateOfCharge else esM.hoursPerTimeStep

        def op2(pyM, loc, compName, p, t):
            rate = getattr(compDict[compName], opRateName)
            return opVar[loc, compName, p, t] == capVar[loc, compName] * rate[loc][p, t] * factor
        setattr(pyM, constrName + '2_' + abbrvName, pyomo.Constraint(constrSet2, pyM.timeSet, rule=op2)) 

def yearlyLimitationConstraint(self, pyM, esM):
        """
        Limit annual commodity imports/exports over the energySystemModel's boundaries for one or multiple
        Source/Sink components.

        :param esM: EnergySystemModel instance representing the energy system in which the component should be modeled.
        :type esM: esM - EnergySystemModel class instance

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        opVar = getattr(pyM, 'op_' + abbrvName)
        limitDict = getattr(pyM, 'yearlyCommodityLimitationDict_' + abbrvName)

        def yearlyLimitationConstraint(pyM, key):
            sumEx = -sum(opVar[loc, compName, p, t] * compDict[compName].sign *
                         esM.periodOccurrences[p]/esM.numberOfYears
                         for loc, compName, p, t in opVar if compName in limitDict[key][1])
            sign = limitDict[key][0]/abs(limitDict[key][0]) if limitDict[key][0] != 0 else 1
            return sign * sumEx <= sign * limitDict[key][0]
        setattr(pyM, 'ConstrYearlyLimitation_' + abbrvName,
                pyomo.Constraint(limitDict.keys(), rule=yearlyLimitationConstraint)) 

def powerFlowDC(self, pyM):
        """
        Ensure that the flow between two locations is equal to the difference between the phase angle variables at
        these locations divided by the reactance of the line between these locations.

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo Concrete Model
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        phaseAngleVar = getattr(pyM, 'phaseAngle_' + self.abbrvName)
        opVar, opVarSet = getattr(pyM, 'op_' + abbrvName), getattr(pyM, 'operationVarSet_' + abbrvName)

        def powerFlowDC(pyM, loc, compName, p, t):
            node1, node2 = compDict[compName]._mapC[loc]
            return (opVar[loc, compName, p, t] - opVar[compDict[compName]._mapI[loc], compName, p, t] ==
                    (phaseAngleVar[node1, compName, p, t]-phaseAngleVar[node2, compName, p, t])/
                    compDict[compName].reactances[loc])
        setattr(pyM, 'ConstrpowerFlowDC_' + abbrvName,  pyomo.Constraint(opVarSet, pyM.timeSet, rule=powerFlowDC)) 

def partLoadOperationOutput(self, pyM):
        """
        Set the required input of a conversion process dependent on the part load efficency.

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo Concrete Model
        """

        compDict, abbrvName = self.componentsDict, self.abbrvName
        discretizationPointConVar = getattr(pyM, 'discretizationPoint_' + self.abbrvName)
        opVar, opVarSet = getattr(pyM, 'op_' + abbrvName), getattr(pyM, 'operationVarSet_' + abbrvName)

        def partLoadOperationOutput(pyM, loc, compName, p, t):        
            nPoints = compDict[compName].nSegments+1
            ### TODO Store the part load levels seperately and do not use 
            # print(list(compDict[compName].discretizedPartLoad.keys()))
            return opVar[loc, compName, p, t] == sum(discretizationPointConVar[loc, compName, discretStep, p, t] * \
                                                 compDict[compName].discretizedPartLoad[list(compDict[compName].discretizedPartLoad.keys())[0]]['xSegments'][discretStep] \
                                                 for discretStep in range(nPoints))
        setattr(pyM, 'ConstrpartLoadOperationOutput_' + abbrvName,  pyomo.Constraint(opVarSet, pyM.timeSet, rule=partLoadOperationOutput)) 

def pointSOS2(self, pyM):
        """
        Ensure that only two consecutive point variables are non-zero while all other point variables are fixed to zero.

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo Concrete Model
        """

        compDict, abbrvName = self.componentsDict, self.abbrvName
        discretizationPointConVar = getattr(pyM, 'discretizationPoint_' + self.abbrvName)
        discretizationSegmentConVar = getattr(pyM, 'discretizationSegmentCon_' + self.abbrvName)
        discretizationPointVarSet = getattr(pyM, 'discretizationPointVarSet_' + self.abbrvName)

        def pointSOS2(pyM, loc, compName, discretStep, p, t):
            points = list(range(compDict[compName].nSegments+1))
            segments = list(range(compDict[compName].nSegments))

            if discretStep == points[0]:
                return discretizationPointConVar[loc, compName, points[0], p, t] <= discretizationSegmentConVar[loc, compName, segments[0], p, t]
            elif discretStep == points[-1]:
                return discretizationPointConVar[loc, compName, points[-1], p, t] <= discretizationSegmentConVar[loc, compName, segments[-1], p, t]
            else:
                return discretizationPointConVar[loc, compName, discretStep, p, t] <= discretizationSegmentConVar[loc, compName, discretStep-1, p, t] + discretizationSegmentConVar[loc, compName, discretStep, p, t]

        setattr(pyM, 'ConstrPointSOS2_' + abbrvName,  pyomo.Constraint(discretizationPointVarSet, pyM.timeSet, rule=pointSOS2)) 

def segmentCapacityConstraint(self, pyM, esM):
        """
        Ensure that the continuous segment variables are in sum equal to the installed capacity of the component.

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo Concrete Model
        """

        compDict, abbrvName = self.componentsDict, self.abbrvName
        discretizationSegmentConVar = getattr(pyM, 'discretizationSegmentCon_' + self.abbrvName)
        capVar = getattr(pyM, 'cap_' + abbrvName)
        opVarSet = getattr(pyM, 'operationVarSet_' + abbrvName)

        def segmentCapacityConstraint(pyM, loc, compName, p, t):
            return sum(discretizationSegmentConVar[loc, compName, discretStep, p, t] for discretStep in range(compDict[compName].nSegments)) == esM.hoursPerTimeStep * capVar[loc, compName]
        setattr(pyM, 'ConstrSegmentCapacity_' + abbrvName,  pyomo.Constraint(opVarSet, pyM.timeSet, rule=segmentCapacityConstraint)) 

def segmentBigM(self, pyM):
        """
        Ensure that the continuous segment variables are zero if the respective binary variable is zero and unlimited otherwise.

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo Concrete Model
        """

        compDict, abbrvName = self.componentsDict, self.abbrvName
        discretizationSegmentConVar = getattr(pyM, 'discretizationSegmentCon_' + self.abbrvName)
        discretizationSegmentBinVar = getattr(pyM, 'discretizationSegmentBin_' + self.abbrvName)
        discretizationSegmentVarSet = getattr(pyM, 'discretizationSegmentVarSet_' + self.abbrvName)

        def segmentBigM(pyM, loc, compName, discretStep, p, t):
            return discretizationSegmentConVar[loc, compName, discretStep, p, t] <= discretizationSegmentBinVar[loc, compName, discretStep, p, t] * compDict[compName].bigM
        setattr(pyM, 'ConstrSegmentBigM_' + abbrvName,  pyomo.Constraint(discretizationSegmentVarSet, pyM.timeSet, rule=segmentBigM)) 

def designBinFix(self, pyM):
        """ 
        Set, if applicable, the installed capacities of a component. 
        
        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel
        """
        compDict, abbrvName, dim = self.componentsDict, self.abbrvName, self.dimension
        designBinVar = getattr(pyM, 'designBin_' + abbrvName)
        designBinVarSet = getattr(pyM, 'designDecisionVarSet_' + abbrvName)

        def designBinFix(pyM, loc, compName):
            return (designBinVar[loc, compName] == compDict[compName].isBuiltFix[loc]
                    if compDict[compName].isBuiltFix is not None else pyomo.Constraint.Skip)
        setattr(pyM, 'ConstrDesignBinFix_' + abbrvName, pyomo.Constraint(designBinVarSet, rule=designBinFix))

    ####################################################################################################################
    #                               Functions for declaring time dependent constraints                                 #
    #################################################################################################################### 

def operationMode1(self, pyM, esM, constrName, constrSetName, opVarName, factorName=None, isStateOfCharge=False):
        """
        Define operation mode 1. The operation [commodityUnit*h] is limited by the installed capacity in:\n
        * [commodityUnit*h] (for storages) or in
        * [commodityUnit] multiplied by the hours per time step (else).\n
        An additional factor can limited the operation further.
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        opVar, capVar = getattr(pyM, opVarName + '_' + abbrvName), getattr(pyM, 'cap_' + abbrvName)
        constrSet1 = getattr(pyM, constrSetName + '1_' + abbrvName)
        factor1 = 1 if isStateOfCharge else esM.hoursPerTimeStep

        def op1(pyM, loc, compName, p, t):
            factor2 = 1 if factorName is None else getattr(compDict[compName], factorName)
            return opVar[loc, compName, p, t] <= factor1 * factor2 * capVar[loc, compName]
        setattr(pyM, constrName + '1_' + abbrvName, pyomo.Constraint(constrSet1, pyM.timeSet, rule=op1)) 

def rampDownMax(self, pyM, esM):
            """
            Ensure that conversion unit is not ramping down too fast by implementing a maximum ramping rate as share of the installed capacity.
            
    
            :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
            :type pyM: pyomo Concrete Model
            """
            compDict, abbrvName = self.componentsDict, self.abbrvName
            
            opVar= getattr(pyM, 'op_' + abbrvName)
            capVar= getattr(pyM, 'cap_' + abbrvName)
            
            constrSetRampDownMax = getattr(pyM,'opConstrSet' + 'rampDownMax_' + abbrvName)
            numberOfTimeSteps = esM.numberOfTimeSteps
    
            def rampDownMax(pyM, loc, compName, p, t):
                rampRateMax = getattr(compDict[compName], 'rampDownMax')
                if (t>=1): # avoid to set constraints twice
                    return (opVar[loc, compName, p, t-1]-opVar[loc, compName, p, t] <= rampRateMax*capVar[loc, compName])
                else:
                    return (opVar[loc, compName, p, numberOfTimeSteps-1]-opVar[loc, compName, p, t] <= rampRateMax*capVar[loc, compName])
            setattr(pyM, 'ConstrRampDownMax_' + abbrvName, pyomo.Constraint(constrSetRampDownMax, pyM.timeSet, rule=rampDownMax)) 

def rampUpMax(self, pyM, esM):
            """
            Ensure that conversion unit is not ramping up too fast by implementing a maximum ramping rate as share of the installed capacity.
            

            :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
            :type pyM: pyomo Concrete Model
            """
            compDict, abbrvName = self.componentsDict, self.abbrvName
            
            opVar= getattr(pyM, 'op_' + abbrvName)
            capVar= getattr(pyM, 'cap_' + abbrvName)
            
            constrSetRampUpMax = getattr(pyM,'opConstrSet' + 'rampUpMax_' + abbrvName)
            numberOfTimeSteps = esM.numberOfTimeSteps
    
            def rampUpMax(pyM, loc, compName, p, t):
                rampRateMax = getattr(compDict[compName], 'rampUpMax')
                if (t>=1): # avoid to set constraints twice
                    return (opVar[loc, compName, p, t]-opVar[loc, compName, p, t-1] <= rampRateMax*capVar[loc, compName])
                else:
                    return (opVar[loc, compName, p, t]-opVar[loc, compName, p, numberOfTimeSteps-1] <= rampRateMax*capVar[loc, compName])
            setattr(pyM, 'ConstrRampUpMax_' + abbrvName, pyomo.Constraint(constrSetRampUpMax, pyM.timeSet, rule=rampUpMax)) 

def operationMode3(self, pyM, esM, constrName, constrSetName, opVarName, opRateName='operationRateMax',
                       isStateOfCharge=False):
        """
        Define operation mode 3. The operation [commodityUnit*h] is limited by an installed capacity multiplied
        with a time series in:\n
        * [commodityUnit*h] (for storages) or in
        * [commodityUnit] multiplied by the hours per time step (else).\n
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        opVar, capVar = getattr(pyM, opVarName + '_' + abbrvName), getattr(pyM, 'cap_' + abbrvName)
        constrSet3 = getattr(pyM, constrSetName + '3_' + abbrvName)
        factor = 1 if isStateOfCharge else esM.hoursPerTimeStep

        def op3(pyM, loc, compName, p, t):
            rate = getattr(compDict[compName], opRateName)
            return opVar[loc, compName, p, t] <= capVar[loc, compName] * rate[loc][p, t] * factor
        setattr(pyM, constrName + '3_' + abbrvName, pyomo.Constraint(constrSet3, pyM.timeSet, rule=op3)) 

def connectSOCs(self, pyM, esM):
        """
        Declare the constraint for connecting the state of charge with the charge and discharge operation:
        the change in the state of charge between two points in time has to match the values of charging and
        discharging (considering the efficiencies of these processes) within the time step in between minus
        the self-discharge of the storage.

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel

        :param esM: EnergySystemModel instance representing the energy system in which the component should be modeled.
        :type esM: esM - EnergySystemModel class instance
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        SOC = getattr(pyM, 'stateOfCharge_' + abbrvName)
        chargeOp, dischargeOp = getattr(pyM, 'chargeOp_' + abbrvName), getattr(pyM, 'dischargeOp_' + abbrvName)
        opVarSet = getattr(pyM, 'operationVarSet_' + abbrvName)

        def connectSOCs(pyM, loc, compName, p, t):
            return (SOC[loc, compName, p, t+1] - SOC[loc, compName, p, t] *
                    (1 - compDict[compName].selfDischarge) ** esM.hoursPerTimeStep ==
                    chargeOp[loc, compName, p, t] * compDict[compName].chargeEfficiency -
                    dischargeOp[loc, compName, p, t] / compDict[compName].dischargeEfficiency)
        setattr(pyM, 'ConstrConnectSOC_' + abbrvName, pyomo.Constraint(opVarSet, pyM.timeSet, rule=connectSOCs)) 

def cyclicLifetime(self, pyM, esM):
        """
        Declare the constraint for limiting the number of full cycle equivalents to stay below cyclic lifetime.

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel

        :param esM: EnergySystemModel instance representing the energy system in which the component should be modeled.
        :type esM: esM - EnergySystemModel class instance
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        chargeOp, capVar = getattr(pyM, 'chargeOp_' + abbrvName), getattr(pyM, 'cap_' + abbrvName)
        capVarSet = getattr(pyM, 'designDimensionVarSet_' + abbrvName)

        def cyclicLifetime(pyM, loc, compName):
            return (sum(chargeOp[loc, compName, p, t] * esM.periodOccurrences[p] for p, t in pyM.timeSet) /
                    esM.numberOfYears <= capVar[loc, compName] *
                    (compDict[compName].stateOfChargeMax - compDict[compName].stateOfChargeMin) *
                    compDict[compName].cyclicLifetime / compDict[compName].economicLifetime[loc]
                    if compDict[compName].cyclicLifetime is not None else pyomo.Constraint.Skip)
        setattr(pyM, 'ConstrCyclicLifetime_' + abbrvName, pyomo.Constraint(capVarSet, rule=cyclicLifetime)) 

def operationMode1_2dim(self, pyM, esM, constrName, constrSetName, opVarName):
        """
        Declare the constraint that the operation [commodityUnit*hour] is limited by the installed
        capacity [commodityUnit] multiplied by the hours per time step.
        Since the flow should either go in one direction or the other, the limitation can be enforced on the sum
        of the forward and backward flow over the line. This leads to one of the flow variables being set to zero
        if a basic solution is obtained during optimization.

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel

        :param esM: EnergySystemModel instance representing the energy system in which the component should be modeled.
        :type esM: esM - EnergySystemModel class instance
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        opVar, capVar = getattr(pyM, opVarName + '_' + abbrvName), getattr(pyM, 'cap_' + abbrvName)
        constrSet1 = getattr(pyM, constrSetName + '1_' + abbrvName)

        def op1(pyM, loc, compName, p, t):
            return opVar[loc, compName, p, t] + opVar[compDict[compName]._mapI[loc], compName, p, t] <= \
                   capVar[loc, compName] * esM.hoursPerTimeStep
        setattr(pyM, constrName + '_' + abbrvName, pyomo.Constraint(constrSet1, pyM.timeSet, rule=op1)) 

def extra_functionality(network,snapshots):

        #Guarantees heat output and electric output nominal powers are proportional
        network.model.chp_nom = Constraint(rule=lambda model : network.links.at["generator","efficiency"]*nom_r*model.link_p_nom["generator"] == network.links.at["boiler","efficiency"]*model.link_p_nom["boiler"])

        #Guarantees c_m p_b1  \leq p_g1
        def backpressure(model,snapshot):
            return c_m*network.links.at["boiler","efficiency"]*model.link_p["boiler",snapshot] <= network.links.at["generator","efficiency"]*model.link_p["generator",snapshot] 
        
        network.model.backpressure = Constraint(list(snapshots),rule=backpressure)
        
        #Guarantees p_g1 +c_v p_b1 \leq p_g1_nom
        def top_iso_fuel_line(model,snapshot):
            return model.link_p["boiler",snapshot] + model.link_p["generator",snapshot] <= model.link_p_nom["generator"]
        
        network.model.top_iso_fuel_line = Constraint(list(snapshots),rule=top_iso_fuel_line) 

def operationModeSOC(self, pyM, esM):
        """
        Declare the constraint that the state of charge [commodityUnit*h] is limited by the installed capacity
        [commodityUnit*h] and the relative maximum state of charge [-].

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel

        :param esM: EnergySystemModel instance representing the energy system in which the component should be modeled.
        :type esM: esM - EnergySystemModel class instance
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        opVar, capVar = getattr(pyM, 'stateOfCharge_' + abbrvName), getattr(pyM, 'cap_' + abbrvName)
        constrSet = getattr(pyM, 'designDimensionVarSet_' + abbrvName)

        # Operation [commodityUnit*h] limited by the installed capacity [commodityUnit*h] multiplied by the relative
        # maximum state of charge.
        def op(pyM, loc, compName, p, t):
            return (opVar[loc, compName, p, t] <=
                    compDict[compName].stateOfChargeMax * capVar[loc, compName])
        setattr(pyM, 'ConstrSOCMaxPrecise_' + abbrvName, pyomo.Constraint(constrSet, pyM.timeSet, rule=op)) 

def intraSOCstart(self, pyM, esM):
        """
        Declare the constraint that the (virtual) state of charge at the beginning of a typical period is zero.

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel

        :param esM: EnergySystemModel instance representing the energy system in which the component should be modeled.
        :type esM: esM - EnergySystemModel class instance
        """
        abbrvName = self.abbrvName
        opVarSet = getattr(pyM, 'operationVarSet_' + abbrvName)
        SOC = getattr(pyM, 'stateOfCharge_' + abbrvName)

        def intraSOCstart(pyM, loc, compName, p):
            return SOC[loc, compName, p, 0] == 0
        setattr(pyM, 'ConstrSOCPeriodStart_' + abbrvName,
                pyomo.Constraint(opVarSet, esM.typicalPeriods, rule=intraSOCstart)) 

def additionalMinPartLoad(self, pyM, esM, constrName, constrSetName, opVarName, opVarBinName, capVarName):
        """
        Set, if applicable, the minimal part load of a component.

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName

        opVar = getattr(pyM, opVarName + '_' + abbrvName)
        opVarBin = getattr(pyM, opVarBinName + '_' + abbrvName)
        capVar = getattr(pyM, capVarName + '_' + abbrvName)
        constrSetMinPartLoad = getattr(pyM, constrSetName + 'partLoadMin_' + abbrvName)
        
        def opMinPartLoad1(pyM, loc, compName, p, t):
            bigM = getattr(compDict[compName], 'bigM')
            return opVar[loc, compName, p, t] <= opVarBin[loc, compName, p, t]*bigM
        setattr(pyM, constrName + 'partLoadMin_1_' + abbrvName, pyomo.Constraint(constrSetMinPartLoad, pyM.timeSet, rule=opMinPartLoad1))
        
        def opMinPartLoad2(pyM, loc, compName, p, t):
            partLoadMin = getattr(compDict[compName], 'partLoadMin')
            bigM = getattr(compDict[compName], 'bigM')
            return opVar[loc, compName, p, t] >= partLoadMin*capVar[loc, compName]-(1-opVarBin[loc, compName, p, t])*bigM
        setattr(pyM, constrName + 'partLoadMin_2_' + abbrvName, pyomo.Constraint(constrSetMinPartLoad, pyM.timeSet, rule=opMinPartLoad2)) 

def bigM(self, pyM):
        """ 
        Enforce the consideration of the binary design variables of a component. 
        
        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        capVar, designBinVar = getattr(pyM, 'cap_' + abbrvName), getattr(pyM, 'designBin_' + abbrvName)
        designBinVarSet = getattr(pyM, 'designDecisionVarSet_' + abbrvName)

        def bigM(pyM, loc, compName):
            return capVar[loc, compName] <= designBinVar[loc, compName] * compDict[compName].bigM
        setattr(pyM, 'ConstrBigM_' + abbrvName, pyomo.Constraint(designBinVarSet, rule=bigM)) 

def capToNbReal(self, pyM):
        """ 
        Determine the components' capacities from the number of installed units.
        
        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel 
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        capVar, nbRealVar = getattr(pyM, 'cap_' + abbrvName), getattr(pyM, 'nbReal_' + abbrvName)
        nbRealVarSet = getattr(pyM, 'continuousDesignDimensionVarSet_' + abbrvName)

        def capToNbReal(pyM, loc, compName):
            return capVar[loc, compName] == nbRealVar[loc, compName] * compDict[compName].capacityPerPlantUnit
        setattr(pyM, 'ConstrCapToNbReal_' + abbrvName, pyomo.Constraint(nbRealVarSet, rule=capToNbReal)) 

def test_instance_constraints(model):
        instance = model.create_instance(report_timing=False)
        for c in instance.component_objects(Constraint):
            c.activate()
            solver = SolverFactory(cfg.solver_name)
            solution = solver.solve(instance)
            if solution.solver.termination_condition == TerminationCondition.infeasible:
                pass
            else:
                print c.name
                c.activate() 

def empty_model(model):
    logger.debug("Storing pyomo model to disk")
    rules = {}
    for obj in model.component_objects(ctype=Constraint):
        if obj.rule is not None:
            rules[obj.name] = obj.rule
            obj.rule = None

    bounds = {}
    for obj in model.component_objects(ctype=Var):
        if obj._bounds_init_rule is not None:
            bounds[obj.name] = obj._bounds_init_rule
            obj._bounds_init_rule = None

    fd, fn = tempfile.mkstemp()
    with os.fdopen(fd, 'wb') as f:
        pickle.dump(model.__getstate__(), f, -1)

    model.__dict__.clear()
    logger.debug("Stored pyomo model to disk")

    gc.collect()
    yield

    logger.debug("Reloading pyomo model")
    with open(fn, 'rb') as f:
        state = pickle.load(f)
    os.remove(fn)
    model.__setstate__(state)

    for n, rule in iteritems(rules):
        getattr(model, n).rule = rule

    for n, bound in iteritems(bounds):
        getattr(model, n)._bounds_init_rule = bound

    logger.debug("Reloaded pyomo model") 

def free_pyomo_initializers(obj):
    obj.construct()
    if isinstance(obj, Var):
        attrs = ('_bounds_init_rule', '_bounds_init_value',
                 '_domain_init_rule', '_domain_init_value',
                 '_value_init_rule', '_value_init_value')
    elif isinstance(obj, Constraint):
        attrs = ('rule', '_init_expr')
    else:
        raise NotImplementedError

    for attr in attrs:
        if hasattr(obj, attr):
            setattr(obj, attr, None) 

def capToNbInt(self, pyM):
        """ 
        Determine the components' capacities from the number of installed units. 
        
        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        capVar, nbIntVar = getattr(pyM, 'cap_' + abbrvName), getattr(pyM, 'nbInt_' + abbrvName)
        nbIntVarSet = getattr(pyM, 'discreteDesignDimensionVarSet_' + abbrvName)

        def capToNbInt(pyM, loc, compName):
            return capVar[loc, compName] == nbIntVar[loc, compName] * compDict[compName].capacityPerPlantUnit
        setattr(pyM, 'ConstrCapToNbInt_' + abbrvName, pyomo.Constraint(nbIntVarSet, rule=capToNbInt)) 

def capacityMinDec(self, pyM):
        """ 
        Enforce the consideration of minimum capacities for components with design decision variables. 
        
        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel
        """
        compDict, abbrvName, dim = self.componentsDict, self.abbrvName, self.dimension
        capVar, designBinVar = getattr(pyM, 'cap_' + abbrvName), getattr(pyM, 'designBin_' + abbrvName)
        designBinVarSet = getattr(pyM, 'designDecisionVarSet_' + abbrvName)

        def capacityMinDec(pyM, loc, compName):
            return (capVar[loc, compName] >= compDict[compName].capacityMin[loc] * designBinVar[loc, compName]
                    if compDict[compName].capacityMin is not None else pyomo.Constraint.Skip)
        setattr(pyM, 'ConstrCapacityMinDec_' + abbrvName, pyomo.Constraint(designBinVarSet, rule=capacityMinDec)) 

def yearlyFullLoadHoursMax(self, pyM, esM):
        """
        Limit the annual full load hours to a maximum value.

        :param esM: EnergySystemModel instance representing the energy system in which the component should be modeled.
        :type esM: esM - EnergySystemModel class instance

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        opVar = getattr(pyM, 'op_' + abbrvName)
        capVar = getattr(pyM, 'cap_' + abbrvName)
        yearlyFullLoadHoursMaxSet = getattr(pyM, 'yearlyFullLoadHoursMaxSet_' + abbrvName)

        def yearlyFullLoadHoursMaxConstraint(pyM, loc, compName, p, t):
            full_load_hours = sum(
                opVar[loc, compName, p, t] * esM.periodOccurrences[p] / esM.numberOfYears for loc, compName, p, t,
                in opVar)
            return full_load_hours <= capVar[loc, compName] * compDict[compName].yearlyFullLoadHoursMax[loc]

        setattr(pyM, 'ConstrYearlyFullLoadHoursMax_' + abbrvName,
                pyomo.Constraint(yearlyFullLoadHoursMaxSet, pyM.timeSet, rule=yearlyFullLoadHoursMaxConstraint))


    ####################################################################################################################
    #  Functions for declaring component contributions to basic energy system constraints and the objective function   #
    #################################################################################################################### 

def cyclicState(self, pyM, esM):
        """
        Declare the constraint for connecting the states of charge: the state of charge at the beginning of a period
        has to be the same as the state of charge in the end of that period.

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel

        :param esM: EnergySystemModel instance representing the energy system in which the component should be modeled.
        :type esM: esM - EnergySystemModel class instance
        """
        abbrvName = self.abbrvName
        opVarSet = getattr(pyM, 'operationVarSet_' + abbrvName)
        SOC = getattr(pyM, 'stateOfCharge_' + abbrvName)
        offsetUp = getattr(pyM, 'stateOfChargeOffsetUp_' + abbrvName)
        offsetDown = getattr(pyM, 'stateOfChargeOffsetDown_' + abbrvName)

        if not pyM.hasTSA:
            def cyclicState(pyM, loc, compName):
                offsetUp_ = offsetUp[loc, compName, 0] if (loc, compName, 0) in offsetUp else 0
                offsetDown_ = offsetDown[loc, compName, 0] if (loc, compName, 0) in offsetDown else 0
                return SOC[loc, compName, 0, 0] == \
                    SOC[loc, compName, 0, esM.timeStepsPerPeriod[-1] + 1] + (offsetUp_ - offsetDown_)
        else:
            SOCInter = getattr(pyM, 'stateOfChargeInterPeriods_' + abbrvName)
            def cyclicState(pyM, loc, compName):
                tLast = esM.interPeriodTimeSteps[-1]
                offsetUp_ = offsetUp[loc, compName, tLast] if (loc, compName, tLast) in offsetUp else 0
                offsetDown_ = offsetDown[loc, compName, tLast] if (loc, compName, tLast) in offsetDown else 0
                return SOCInter[loc, compName, 0] == \
                    SOCInter[loc, compName, tLast] + (offsetUp_ - offsetDown_)
        setattr(pyM, 'ConstrCyclicState_' + abbrvName, pyomo.Constraint(opVarSet, rule=cyclicState)) 

def yearlyFullLoadHoursMin(self, pyM, esM):
        # TODO: Add deprecation warning to sourceSink.yearlyLimitConstraint and call this function in it
        # TODO: Use a set to declare Constraint only for components with full load hour limit set.
        """
        Limit the annual full load hours to a minimum value.

        :param esM: EnergySystemModel instance representing the energy system in which the component should be modeled.
        :type esM: esM - EnergySystemModel class instance

        :param pyM: pyomo ConcreteModel which stores the mathematical formulation of the model.
        :type pyM: pyomo ConcreteModel
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        opVar = getattr(pyM, 'op_' + abbrvName)
        capVar = getattr(pyM, 'cap_' + abbrvName)
        yearlyFullLoadHoursMinSet = getattr(pyM, 'yearlyFullLoadHoursMinSet_' + abbrvName)


        def yearlyFullLoadHoursMinConstraint(pyM, loc, compName, p, t):
            full_load_hours = sum(
                opVar[loc, compName, p, t] * esM.periodOccurrences[p] / esM.numberOfYears for loc, compName, p, t,
                in opVar)
            return full_load_hours >= capVar[loc, compName] * compDict[compName].yearlyFullLoadHoursMin[loc]

        setattr(pyM, 'ConstrYearlyFullLoadHoursMin_' + abbrvName,
                pyomo.Constraint(yearlyFullLoadHoursMinSet, pyM.timeSet, rule=yearlyFullLoadHoursMinConstraint)) 

def operationMode4(self, pyM, esM, constrName, constrSetName, opVarName, opRateName='operationRateFix'):
        """
        Define operation mode 4. The operation [commodityUnit*h] is equal to a time series in.
        """
        compDict, abbrvName = self.componentsDict, self.abbrvName
        opVar = getattr(pyM, opVarName + '_' + abbrvName)
        constrSet4 = getattr(pyM, constrSetName + '4_' + abbrvName)

        def op4(pyM, loc, compName, p, t):
            rate = getattr(compDict[compName], opRateName)
            return opVar[loc, compName, p, t] == rate[loc][p, t]
        setattr(pyM, constrName + '4_' + abbrvName, pyomo.Constraint(constrSet4, pyM.timeSet, rule=op4))