def get_capacity_constraint(backend_model, parameter, loc_tech,
                            _equals=None, _max=None, _min=None, scale=None):

    decision_variable = getattr(backend_model, parameter)

    if not _equals:
        _equals = get_param(backend_model, parameter + '_equals', loc_tech)
    if not _max:
        _max = get_param(backend_model, parameter + '_max', loc_tech)
    if not _min:
        _min = get_param(backend_model, parameter + '_min', loc_tech)
    if po.value(_equals) is not False and po.value(_equals) is not None:
        if np.isinf(po.value(_equals)):
            e = exceptions.ModelError
            raise e('Cannot use inf for {}_equals for loc:tech `{}`'.format(parameter, loc_tech))
        if scale:
            _equals *= scale
        return decision_variable[loc_tech] == _equals
    else:
        if po.value(_min) == 0 and np.isinf(po.value(_max)):
            return po.Constraint.NoConstraint
        else:
            if scale:
                _max *= scale
                _min *= scale
            return (_min, decision_variable[loc_tech], _max)

    def solveModel(self, x, y, z):
        model = self.model
        opt = SolverFactory(self.config.solver)
        opt.options.update(self.config.solver_options)

        results = opt.solve(
            model, keepfiles=self.keepfiles, tee=self.stream_solver)

        if ((results.solver.status == SolverStatus.ok)
                and (results.solver.termination_condition == TerminationCondition.optimal)):
            model.solutions.load_from(results)
            for i in range(0, self.lx):
                x[i] = value(self.TRF.xvars[i])
            for i in range(0, self.ly):
                y[i] = value(self.TRF.y[i+1])
            for i in range(0, self.lz):
                z[i] = value(self.TRF.zvars[i])

            for obj in model.component_data_objects(Objective,active=True):
                return True, obj()

        else:
            print("Waring: solver Status: " + str(results.solver.status))
            print("And Termination Conditions: " + str(results.solver.termination_condition))
            return False, 0

def disjunctive_bound(var, scope):
    """Compute the disjunctive bounds for a variable in a given scope.

    Args:
        var (_VarData): Variable for which to compute bound
        scope (Component): The scope in which to compute the bound. If not a
            _DisjunctData, it will walk up the tree and use the scope of the
            most immediate enclosing _DisjunctData.

    Returns:
        numeric: the tighter of either the disjunctive lower bound, the
            variable lower bound, or (-inf, inf) if neither exist.

    """
    # Initialize to the global variable bound
    var_bnd = (
        value(var.lb) if var.has_lb() else -inf,
        value(var.ub) if var.has_ub() else inf)
    possible_disjunct = scope
    while possible_disjunct is not None:
        try:
            disj_bnd = possible_disjunct._disj_var_bounds.get(var, (-inf, inf))
            disj_bnd = (
                max(var_bnd[0], disj_bnd[0]),
                min(var_bnd[1], disj_bnd[1]))
            return disj_bnd
        except AttributeError:
            # possible disjunct does not have attribute '_disj_var_bounds'.
            # Try again with the scope's parent block.
            possible_disjunct = possible_disjunct.parent_block()
    # Unable to find '_disj_var_bounds' attribute within search scope.
    return var_bnd

 def tear_diff_direct(self, G, tears):
     """
     Returns numpy arrays of values for src and dest members
     for all edges in the tears list of edge indexes.
     """
     svals = []
     dvals = []
     edge_list = self.idx_to_edge(G)
     for tear in tears:
         arc = G.edges[edge_list[tear]]["arc"]
         src, dest = arc.src, arc.dest
         sf = arc.expanded_block.component("splitfrac")
         for name, mem in src.iter_vars(names=True):
             if src.is_extensive(name) and sf is not None:
                 # TODO: same as above, what if there's no splitfrac
                 svals.append(value(mem * sf))
             else:
                 svals.append(value(mem))
             try:
                 index = mem.index()
             except AttributeError:
                 index = None
             dvals.append(value(self.source_dest_peer(arc, name, index)))
     svals = numpy.array(svals)
     dvals = numpy.array(dvals)
     return svals, dvals

def make2dPlot(expr, numticks=10, show_plot=False):
    mc_ccVals = [None] * (numticks + 1)
    mc_cvVals = [None] * (numticks + 1)
    aff_cc = [None] * (numticks + 1)
    aff_cv = [None] * (numticks + 1)
    fvals = [None] * (numticks + 1)
    mc_expr = mc(expr)
    x = next(identify_variables(expr))  # get the first variable
    tick_length = (x.ub - x.lb) / numticks
    xaxis = [x.lb + tick_length * n for n in range(numticks + 1)]

    x_val = value(x)  # initial value of x
    cc = mc_expr.subcc()  # Concave overestimator subgradient at x_val
    cv = mc_expr.subcv()  # Convex underestimator subgradient at x_val
    f_cc = mc_expr.concave()  # Concave overestimator value at x_val
    f_cv = mc_expr.convex()  # Convex underestimator value at x_val
    for i, x_tick in enumerate(xaxis):
        aff_cc[i] = cc[x] * (x_tick - x_val) + f_cc
        aff_cv[i] = cv[x] * (x_tick - x_val) + f_cv
        mc_expr.changePoint(x, x_tick)
        mc_ccVals[i] = mc_expr.concave()
        mc_cvVals[i] = mc_expr.convex()
        fvals[i] = value(expr)
    if show_plot:
        import matplotlib.pyplot as plt
        plt.plot(xaxis, fvals, 'r', xaxis, mc_ccVals, 'b--', xaxis,
                 mc_cvVals, 'b--', xaxis, aff_cc, 'k|', xaxis, aff_cv, 'k|')
        plt.show()
    return mc_ccVals, mc_cvVals, aff_cc, aff_cv

def copy_var_list_values(from_list, to_list, config, skip_stale=False):
    """Copy variable values from one list to another."""
    for v_from, v_to in zip(from_list, to_list):
        if skip_stale and v_from.stale:
            continue  # Skip stale variable values.
        try:
            v_to.set_value(value(v_from, exception=False))
            if skip_stale:
                v_to.stale = False
        except ValueError as err:
            err_msg = getattr(err, 'message', str(err))
            var_val = value(v_from)
            rounded_val = int(round(var_val))
            # Check to see if this is just a tolerance issue
            if 'is not in domain Binary' in err_msg and (
                    fabs(var_val - 1) <= config.integer_tolerance or
                    fabs(var_val) <= config.integer_tolerance):
                v_to.set_value(rounded_val)
            elif 'is not in domain Integers' in err_msg and (
                    fabs(var_val - rounded_val) <= config.integer_tolerance):
                v_to.set_value(rounded_val)
            # Value is zero, but shows up as slightly less than zero.
            elif 'is not in domain NonNegativeReals' in err_msg and (
                    fabs(var_val) <= config.zero_tolerance):
                v_to.set_value(0)
            else:
                raise

    def _estimate_M(self, expr, name):
        # Calculate a best guess at M
        repn = generate_standard_repn(expr)
        M = [0, 0]

        if not repn.is_nonlinear():
            if repn.constant is not None:
                for i in (0, 1):
                    if M[i] is not None:
                        M[i] += repn.constant

            for i, coef in enumerate(repn.linear_coefs or []):
                var = repn.linear_vars[i]
                bounds = (value(var.lb), value(var.ub))
                for i in (0, 1):
                    # reverse the bounds if the coefficient is negative
                    if coef > 0:
                        j = i
                    else:
                        j = 1 - i

                    if bounds[i] is not None:
                        M[j] += value(bounds[i]) * coef
                    else:
                        raise GDP_Error(
                            "Cannot estimate M for "
                            "expressions with unbounded variables."
                            "\n\t(found unbounded var %s while processing "
                            "constraint %s)" % (var.name, name))
        else:
            raise GDP_Error("Cannot estimate M for nonlinear "
                            "expressions.\n\t(found while processing "
                            "constraint %s)" % name)

        return tuple(M)

def unit_capacity_systemwide_constraint_rule(backend_model, tech):
    """
    Set constraints to limit the number of purchased units of a single technology
    type across all locations in the model.

    The first valid case is applied:

    .. container:: scrolling-wrapper

        .. math::

            \\sum_{loc}\\boldsymbol{units}(loc::tech) + \\boldsymbol{purchased}(loc::tech)
            \\begin{cases}
                = units_{equals, systemwide}(tech),&
                    \\text{if } units_{equals, systemwide}(tech)\\\\
                \\leq units_{max, systemwide}(tech),&
                    \\text{if } units_{max, systemwide}(tech)\\\\
                \\text{unconstrained},& \\text{otherwise}
            \\end{cases}
            \\forall tech \\in techs

    """

    if tech in backend_model.techs_transmission_names:
        all_loc_techs = [
            i for i in backend_model.loc_techs_transmission
            if i.split('::')[1].split(':')[0] == tech
        ]
        multiplier = 2  # there are always two technologies associated with one link
    else:
        all_loc_techs = [
            i for i in backend_model.loc_techs
            if i.split('::')[1] == tech
        ]
        multiplier = 1

    max_systemwide = get_param(backend_model, 'units_max_systemwide', tech)
    equals_systemwide = get_param(backend_model, 'units_equals_systemwide', tech)

    if np.isinf(po.value(max_systemwide)) and not equals_systemwide:
        return po.Constraint.NoConstraint
    elif equals_systemwide and np.isinf(po.value(equals_systemwide)):
        raise ValueError(
            'Cannot use inf for energy_cap_equals_systemwide for tech `{}`'.format(tech)
        )

    sum_expr_units = sum(
        backend_model.units[loc_tech] for loc_tech in all_loc_techs
        if loc_tech_is_in(backend_model, loc_tech, 'loc_techs_milp')
    )
    sum_expr_purchase = sum(
        backend_model.purchased[loc_tech] for loc_tech in all_loc_techs
        if loc_tech_is_in(backend_model, loc_tech, 'loc_techs_purchase')
    )

    if equals_systemwide:
        return sum_expr_units + sum_expr_purchase == equals_systemwide * multiplier
    else:
        return sum_expr_units + sum_expr_purchase <= max_systemwide * multiplier

def resource_availability_supply_plus_constraint_rule(backend_model, loc_tech, timestep):
    """
    Limit production from supply_plus techs to their available resource.

    .. container:: scrolling-wrapper

        .. math::

            \\boldsymbol{resource_{con}}(loc::tech, timestep)
            \\leq available\\_resource(loc::tech, timestep)
            \\quad \\forall loc::tech \\in loc::techs_{supply^{+}}, \\forall timestep \\in timesteps

    If :math:`force\\_resource(loc::tech)` is set:

    .. container:: scrolling-wrapper

        .. math::

            \\boldsymbol{resource_{con}}(loc::tech, timestep)
            = available\\_resource(loc::tech, timestep)
            \\quad \\forall loc::tech \\in loc::techs_{supply^{+}}, \\forall timestep \\in timesteps

    Where:

    .. container:: scrolling-wrapper

        .. math::

            available\\_resource(loc::tech, timestep) = resource(loc::tech, timestep)
            \\times resource_{scale}(loc::tech)

    if :math:`loc::tech` is in :math:`loc::techs_{area}`:

    .. container:: scrolling-wrapper

        .. math::

            available\\_resource(loc::tech, timestep) = resource(loc::tech, timestep)
            \\times resource_{scale}(loc::tech)
            \\times resource_{area}(loc::tech)

    """
    resource = get_param(backend_model, 'resource', (loc_tech, timestep))
    resource_scale = get_param(backend_model, 'resource_scale', loc_tech)
    force_resource = get_param(backend_model, 'force_resource', loc_tech)
    resource_unit = get_param(backend_model, 'resource_unit', loc_tech)

    if po.value(resource_unit) == 'energy_per_area':
        available_resource = resource * resource_scale * backend_model.resource_area[loc_tech]
    elif po.value(resource_unit) == 'energy_per_cap':
        available_resource = resource * resource_scale * backend_model.energy_cap[loc_tech]
    else:
        available_resource = resource * resource_scale

    if po.value(force_resource):
        return backend_model.resource_con[loc_tech, timestep] == available_resource
    else:
        return backend_model.resource_con[loc_tech, timestep] <= available_resource

def log_infeasible_constraints(m, tol=1E-6, logger=logger):
    """Print the infeasible constraints in the model.

    Uses the current model state. Uses pyomo.util.infeasible logger unless one
    is provided.

    Args:
        m (Block): Pyomo block or model to check
        tol (float): feasibility tolerance

    """
    for constr in m.component_data_objects(
            ctype=Constraint, active=True, descend_into=True):
        # if constraint is an equality, handle differently
        if constr.equality and fabs(value(constr.lower - constr.body)) >= tol:
            logger.info('CONSTR {}: {} != {}'.format(
                constr.name, value(constr.body), value(constr.lower)))
            continue
        # otherwise, check LB and UB, if they exist
        if constr.has_lb() and value(constr.lower - constr.body) >= tol:
            logger.info('CONSTR {}: {} < {}'.format(
                constr.name, value(constr.body), value(constr.lower)))
        if constr.has_ub() and value(constr.body - constr.upper) >= tol:
            logger.info('CONSTR {}: {} > {}'.format(
                constr.name, value(constr.body), value(constr.upper)))

 def getInitialValue(self):
     x = np.zeros(self.lx, dtype=float)
     y = np.zeros(self.ly, dtype=float)
     z = np.zeros(self.lz, dtype=float)
     for i in range(0, self.lx):
         x[i] = value(self.TRF.xvars[i])
     for i in range(0, self.ly):
         #initialization of y?
         y[i] = 1
     for i in range(0, self.lz):
         z[i] = value(self.TRF.zvars[i])
     return x, y, z

def energy_capacity_systemwide_constraint_rule(backend_model, tech):
    """
    Set constraints to limit the capacity of a single technology type across all locations in the model.

    The first valid case is applied:

    .. container:: scrolling-wrapper

        .. math::

            \\sum_{loc}\\boldsymbol{energy_{cap}}(loc::tech)
            \\begin{cases}
                = energy_{cap, equals, systemwide}(loc::tech),&
                    \\text{if } energy_{cap, equals, systemwide}(loc::tech)\\\\
                \\leq energy_{cap, max, systemwide}(loc::tech),&
                    \\text{if } energy_{cap, max, systemwide}(loc::tech)\\\\
                \\text{unconstrained},& \\text{otherwise}
            \\end{cases}
            \\forall tech \\in techs

    """

    if tech in backend_model.techs_transmission_names:
        all_loc_techs = [
            i for i in backend_model.loc_techs_transmission
            if i.split('::')[1].split(':')[0] == tech
        ]
        multiplier = 2  # there are always two technologies associated with one link
    else:
        all_loc_techs = [
            i for i in backend_model.loc_techs
            if i.split('::')[1] == tech
        ]
        multiplier = 1

    max_systemwide = get_param(backend_model, 'energy_cap_max_systemwide', tech)
    equals_systemwide = get_param(backend_model, 'energy_cap_equals_systemwide', tech)

    if np.isinf(po.value(max_systemwide)) and not equals_systemwide:
        return po.Constraint.NoConstraint
    elif equals_systemwide and np.isinf(po.value(equals_systemwide)):
        raise exceptions.ModelError(
            'Cannot use inf for energy_cap_equals_systemwide for tech `{}`'.format(tech)
        )

    sum_expr = sum(backend_model.energy_cap[loc_tech] for loc_tech in all_loc_techs)

    if equals_systemwide:
        return sum_expr == equals_systemwide * multiplier
    else:
        return sum_expr <= max_systemwide * multiplier

def solve_NLP_feas(solve_data, config):
    m = solve_data.working_model.clone()
    add_feas_slacks(m)
    MindtPy = m.MindtPy_utils
    next(m.component_data_objects(Objective, active=True)).deactivate()
    for constr in m.component_data_objects(
            ctype=Constraint, active=True, descend_into=True):
        constr.deactivate()
    MindtPy.MindtPy_feas.activate()
    MindtPy.MindtPy_feas_obj = Objective(
        expr=sum(s for s in MindtPy.MindtPy_feas.slack_var[...]),
        sense=minimize)
    for v in MindtPy.variable_list:
        if v.is_binary():
            v.fix(int(round(v.value)))
    # m.pprint()  #print nlp feasibility problem for debugging
    with SuppressInfeasibleWarning():
        feas_soln = SolverFactory(config.nlp_solver).solve(
            m, **config.nlp_solver_args)
    subprob_terminate_cond = feas_soln.solver.termination_condition
    if subprob_terminate_cond is tc.optimal:
        copy_var_list_values(
            MindtPy.variable_list,
            solve_data.working_model.MindtPy_utils.variable_list,
            config)
        pass
    elif subprob_terminate_cond is tc.infeasible:
        raise ValueError('Feasibility NLP infeasible. '
                         'This should never happen.')
    else:
        raise ValueError(
            'MindtPy unable to handle feasibility NLP termination condition '
            'of {}'.format(subprob_terminate_cond))

    var_values = [v.value for v in MindtPy.variable_list]
    duals = [0 for _ in MindtPy.constraint_list]

    for i, constr in enumerate(MindtPy.constraint_list):
        # TODO rhs only works if constr.upper and constr.lower do not both have values.
        # Sometimes you might have 1 <= expr <= 1. This would give an incorrect rhs of 2.
        rhs = ((0 if constr.upper is None else constr.upper) +
               (0 if constr.lower is None else constr.lower))
        sign_adjust = 1 if value(constr.upper) is None else -1
        duals[i] = sign_adjust * max(
            0, sign_adjust * (rhs - value(constr.body)))

    if value(MindtPy.MindtPy_feas_obj.expr) == 0:
        raise ValueError(
            'Problem is not feasible, check NLP solver')

    return var_values, duals

def cost_var_constraint_rule(backend_model, cost, loc_tech, timestep):
    """
    Calculate costs from time-varying decision variables

    .. container:: scrolling-wrapper

        .. math::

            \\boldsymbol{cost_{var}}(cost, loc::tech, timestep) = cost_{prod}(cost, loc::tech, timestep) + cost_{con}(cost, loc::tech, timestep)

            cost_{prod}(cost, loc::tech, timestep) = cost_{om\\_prod}(cost, loc::tech, timestep) \\times weight(timestep) \\times \\boldsymbol{carrier_{prod}}(loc::tech::carrier, timestep)

            prod\\_con\\_eff =
            \\begin{cases}
                = \\boldsymbol{resource_{con}}(loc::tech, timestep),&
                    \\text{if } loc::tech \\in loc\\_techs\\_supply\\_plus \\\\
                = \\frac{\\boldsymbol{carrier_{prod}}(loc::tech::carrier, timestep)}{energy_eff(loc::tech, timestep)},&
                    \\text{if } loc::tech \\in loc\\_techs\\_supply \\\\
            \\end{cases}

            cost_{con}(cost, loc::tech, timestep) = cost_{om\\_con}(cost, loc::tech, timestep) \\times weight(timestep) \\times prod\\_con\\_eff

    """
    model_data_dict = backend_model.__calliope_model_data__

    cost_om_prod = get_param(backend_model, 'cost_om_prod', (cost, loc_tech, timestep))
    cost_om_con = get_param(backend_model, 'cost_om_con', (cost, loc_tech, timestep))
    weight = backend_model.timestep_weights[timestep]

    loc_tech_carrier = model_data_dict['data']['lookup_loc_techs'][loc_tech]

    if po.value(cost_om_prod):
        cost_prod = cost_om_prod * weight * backend_model.carrier_prod[loc_tech_carrier, timestep]
    else:
        cost_prod = 0

    if loc_tech_is_in(backend_model, loc_tech, 'loc_techs_supply_plus') and cost_om_con:
        cost_con = cost_om_con * weight * backend_model.resource_con[loc_tech, timestep]
    elif loc_tech_is_in(backend_model, loc_tech, 'loc_techs_supply') and cost_om_con:
        energy_eff = get_param(backend_model, 'energy_eff', (loc_tech, timestep))
        if po.value(energy_eff) > 0:  # in case energy_eff is zero, to avoid an infinite value
            cost_con = cost_om_con * weight * (backend_model.carrier_prod[loc_tech_carrier, timestep] / energy_eff)
        else:
            cost_con = 0
    else:
        cost_con = 0

    backend_model.cost_var_rhs[cost, loc_tech, timestep].expr = cost_prod + cost_con
    return (backend_model.cost_var[cost, loc_tech, timestep] ==
            backend_model.cost_var_rhs[cost, loc_tech, timestep])

 def generate_gofx(self, G, tears):
     edge_list = self.idx_to_edge(G)
     gofx = []
     for tear in tears:
         arc = G.edges[edge_list[tear]]["arc"]
         src = arc.src
         sf = arc.expanded_block.component("splitfrac")
         for name, mem in src.iter_vars(names=True):
             if src.is_extensive(name) and sf is not None:
                 # TODO: same as above, what if there's no splitfrac
                 gofx.append(value(mem * sf))
             else:
                 gofx.append(value(mem))
     gofx = numpy.array(gofx)
     return gofx

 def pass_single_value(self, port, name, member, val, fixed):
     """
     Fix the value of the port member and add it to the fixed set.
     If the member is an expression, appropriately fix the value of
     its free variable. Error if the member is already fixed but
     different from val, or if the member has more than one free
     variable."
     """
     eq_tol = self.options["almost_equal_tol"]
     if member.is_fixed():
         if abs(value(member) - val) > eq_tol:
             raise RuntimeError(
                 "Member '%s' of port '%s' is already fixed but has a "
                 "different value (by > %s) than what is being passed to it"
                 % (name, port.name, eq_tol))
     elif member.is_expression_type():
         repn = generate_standard_repn(member - val)
         if repn.is_linear() and len(repn.linear_vars) == 1:
             # fix the value of the single variable
             fval = (0 - repn.constant) / repn.linear_coefs[0]
             var = repn.linear_vars[0]
             fixed.add(var)
             var.fix(fval)
         else:
             raise RuntimeError(
                 "Member '%s' of port '%s' had more than "
                 "one free variable when trying to pass a value "
                 "to it. Please fix more variables before passing "
                 "to this port." % (name, port.name))
     else:
         fixed.add(member)
         member.fix(val)

    def subproblem_solve(gdp, config):
        subproblem = gdp.clone()
        TransformationFactory('gdp.bigm').apply_to(subproblem)
        main_obj = next(subproblem.component_data_objects(Objective, active=True))
        obj_sign = 1 if main_obj.sense == minimize else -1

        try:
            result = SolverFactory(config.solver).solve(subproblem, **config.solver_args)
        except RuntimeError as e:
            config.logger.warning(
                "Solver encountered RuntimeError. Treating as infeasible. "
                "Msg: %s\n%s" % (str(e), traceback.format_exc()))
            var_values = [v.value for v in subproblem.GDPbb_utils.variable_list]
            return obj_sign * float('inf'), SolverResults(), var_values

        var_values = [v.value for v in subproblem.GDPbb_utils.variable_list]
        term_cond = result.solver.termination_condition
        if result.solver.status is SolverStatus.ok and any(
                term_cond == valid_cond for valid_cond in (tc.optimal, tc.locallyOptimal, tc.feasible)):
            return value(main_obj.expr), result, var_values
        elif term_cond == tc.unbounded:
            return obj_sign * float('-inf'), result, var_values
        elif term_cond == tc.infeasible:
            return obj_sign * float('inf'), result, var_values
        else:
            config.logger.warning("Unknown termination condition of %s" % term_cond)
            return obj_sign * float('inf'), result, var_values

    def register_var(self, var, lb, ub):
        """Registers a new variable."""
        var_idx = self.var_to_idx[var]
        inf = float('inf')

        # Guard against errant None values in lb and ub
        lb = -inf if lb is None else lb
        ub = inf if ub is None else ub

        lb = max(var.lb if var.has_lb() else -inf, lb)
        ub = min(var.ub if var.has_ub() else inf, ub)
        var_val = value(var, exception=False)
        if lb == -inf:
            lb = -500000
            logger.warning(
                'Var %s missing lower bound. Assuming LB of %s'
                % (var.name, lb))
        if ub == inf:
            ub = 500000
            logger.warning(
                'Var %s missing upper bound. Assuming UB of %s'
                % (var.name, ub))
        if var_val is None:
            var_val = (lb + ub) / 2
            self.missing_value_warnings.append(
                'Var %s missing value. Assuming midpoint value of %s'
                % (var.name, var_val))
        return self.mcpp.newVar(
            lb, var_val, ub, self.num_vars, var_idx)

def energy_capacity_max_purchase_constraint_rule(backend_model, loc_tech):
    """
    Set maximum energy capacity decision variable upper bound as a function of
    binary purchase variable

    The first valid case is applied:

    .. container:: scrolling-wrapper

        .. math::

            \\frac{\\boldsymbol{energy_{cap}}(loc::tech)}{energy_{cap, scale}(loc::tech)}
            \\begin{cases}
                = energy_{cap, equals}(loc::tech) \\times \\boldsymbol{purchased}(loc::tech),&
                    \\text{if } energy_{cap, equals}(loc::tech)\\\\
                \\leq energy_{cap, max}(loc::tech) \\times \\boldsymbol{purchased}(loc::tech),&
                    \\text{if } energy_{cap, max}(loc::tech)\\\\
                \\text{unconstrained},& \\text{otherwise}
            \\end{cases}
            \\forall loc::tech \\in loc::techs_{purchase}

    """
    energy_cap_max = get_param(backend_model, 'energy_cap_max', loc_tech)
    energy_cap_equals = get_param(backend_model, 'energy_cap_equals', loc_tech)
    energy_cap_scale = get_param(backend_model, 'energy_cap_scale', loc_tech)

    if po.value(energy_cap_equals):
        return backend_model.energy_cap[loc_tech] == (
            energy_cap_equals * energy_cap_scale * backend_model.purchased[loc_tech]
        )

    else:
        return backend_model.energy_cap[loc_tech] <= (
            energy_cap_max * energy_cap_scale * backend_model.purchased[loc_tech]
        )

 def _collect_bilinear(self, expr, bilin, quad):
     if not expr.is_expression():
         return
     if type(expr) is _ProductExpression:
         if len(expr._numerator) != 2:
             for e in expr._numerator:
                 self._collect_bilinear(e, bilin, quad)
             # No need to check denominator, as this is poly_degree==2
             return
         if not isinstance(expr._numerator[0], _VarData) or \
                 not isinstance(expr._numerator[1], _VarData):
             raise RuntimeError("Cannot yet handle complex subexpressions")
         if expr._numerator[0] is expr._numerator[1]:
             quad.append( (expr, expr._numerator[0]) )
         else:
             bilin.append( (expr, expr._numerator[0], expr._numerator[1]) )
         return
     if type(expr) is _PowExpression and value(expr._args[1]) == 2:
         # Note: directly testing the value of the exponent above is
         # safe: we have already verified that this expression is
         # polynominal, so the exponent must be constant.
         tmp = _ProductExpression()
         tmp._numerator = [ expr._args[0], expr._args[0] ]
         tmp._denominator = []
         expr._args = (tmp, as_numeric(1))
         #quad.append( (tmp, tmp._args[0]) )
         self._collect_bilinear(tmp, bilin, quad)
         return 
     # All other expression types
     for e in expr._args:
         self._collect_bilinear(e, bilin, quad)