How to use the pyomo.environ function in Pyomo

def test_create_min_rbf_model(self):
        """Test the create_min_rbf_model function.

        This test simply checks whether the function returns a valid
        pyomo.ConcreteModel object.
        """
        model = d1.create_min_rbf_model(self.settings, self.n, self.k,
                                        self.var_lower, self.var_upper,
                                        self.integer_vars, None, self.node_pos,
                                        self.rbf_lambda, self.rbf_h)
        self.assertIsInstance(model, pyomo.environ.ConcreteModel)
        model = d1.create_min_rbf_model(
            self.settings, 10, 20, np.array([0] * 10),np.array([1] * 10),
            np.array([i for i in range(10)]),
            (np.array([0]), np.array([]),
             [(0, 0, np.array([i for i in range(10)]))]),
            np.random.randint(0, 2, size=(20, 10)),
            np.random.uniform(size=20), np.random.uniform(size=11))
        self.assertIsInstance(model, pyomo.environ.ConcreteModel)

def _init(M):
        t0 = M.t.first()
        yield M.x[1, t0] == 0
        yield M.x[2, t0] == 1
        yield aml.ConstraintList.End

def __init__(self, model):
        """

        Parameters
        ----------
        model : ConcreteModel
            Pyomo concrete model
        """
        temporal_dir = tempfile.mkdtemp()
        try:
            filename = os.path.join(temporal_dir, "pynumero_pyomo")
            objectives = model.component_map(aml.Objective, active=True)
            if len(objectives) == 0:
                model._dummy_obj = aml.Objective(expr=0.0)

            model.write(filename+'.nl', 'nl', io_options={"symbolic_solver_labels": True})

            fname, symbolMap = pyomo.opt.WriterFactory('nl')(model, filename, lambda x:True, {})
            varToIndex = pyomo.core.kernel.component_map.ComponentMap()
            conToIndex = pyomo.core.kernel.component_map.ComponentMap()
            for name, obj in six.iteritems(symbolMap.bySymbol):
                if name[0] == 'v':
                    varToIndex[obj()] = int(name[1:])
                elif name[0] == 'c':
                    conToIndex[obj()] = int(name[1:])

            self._varToIndex = varToIndex
            self._conToIndex = conToIndex

# set list of available diameters
    diameters = dic_diam_costs.keys()

    # build concrete pyomo model
    model = py.ConcreteModel()

    # sets for nodes, arcs, diameters, scenarios
    model.nodes = py.Set(initialize=graph.nodes)
    model.arcs = py.Set(initialize=list(distances.keys()), dimen=2)
    # diameters assuming that each pipe has the same diameter options
    model.diameters = py.Set(initialize=diameters)
    # if we have special scenarios, scenario names are tuples, otherwise not
    if robust:
        # set indices for each scenario by its nodePair = (startnode, endnode)
        model.scenarios = py.Set(initialize=specialScenarioNames, dimen=2)
    else:
        # set indices for each timeStep number
        model.scenarios = py.Set(initialize=specialScenarioNames, dimen=1)

    # create variables binaries x are the same for each scenario
    # pressure variables are different for each scenario
    model.x = py.Var(model.arcs, model.diameters, domain=py.Binary)
    if robust:
        def pressureBounds(model, node, startnode, endnode):
            return dic_node_minPress[node] ** 2, dic_node_maxPress[node] ** 2

        model.pi = py.Var(model.nodes, model.scenarios, bounds=pressureBounds)
    else:
        def pressureBounds(model, node, timeStep):
            return dic_node_minPress[node] ** 2, dic_node_maxPress[node] ** 2

def create_relaxation(a, b):
    m = pe.ConcreteModel()
    m.x = pe.Var(bounds=(-20.0, 20.0))
    m.x_sq = pe.Var()
    m.y = pe.Var(bounds=(-20.0, 20.0))
    m.z = pe.Var()

    m.objective = pe.Objective(expr=(a - m.x)**2 + b*m.z**2)
    m.con1 = pe.Constraint(expr=m.z == m.y - m.x_sq)
    m.x_sq_con = coramin.relaxations.PWXSquaredRelaxation()
    m.x_sq_con.build(x=m.x, aux_var=m.x_sq, use_linear_relaxation=True)

    return m

def _get_subs(self, instance):
        subs = {} # Substitute one var for another from a * x + b * y + c = 0
        subs_map = {} # id -> var
        fixes = [] # fix a variable from a * x + c = 0
        cnstr = set() # constraints to deactivate
        rset = set() # variables used in a sub or fixed
        for c in instance.component_data_objects(pyo.Constraint, active=True):
            if (
                pyo.value(c.lower) is not None and
                pyo.value(c.lower) == pyo.value(c.upper) and
                c.body.polynomial_degree() == 1
            ):
                repn = generate_standard_repn(c.body)
                assert len(repn.nonlinear_vars) == 0
                assert len(repn.quadratic_vars) == 0
                if len(repn.linear_vars) > 2:
                    continue
                elif len(repn.linear_vars) < 1:
                    _log.warning("Constraint with no vars {}: {}".format(c, c.expr))
                    continue
                b0 = repn.constant - pyo.value(c.upper)
                v0 = repn.linear_vars[0]
                a0 = repn.linear_coefs[0]

def build(self):
        """
        Add model equations to the unit model.  This is called by a default block
        construnction rule when the unit model is created.
        """
        super().build() # Basic unit model build/read config
        config = self.config # shorter config pointer

        # The thermodynamic expression writer object, te, writes expressions
        # including external function calls to calculate thermodynamic quantities
        # from a set of state variables.
        _assert_properties(config.property_package)
        te = ThermoExpr(blk=self, parameters=config.property_package)

        eff = self.efficiency_isentropic = pyo.Var(
            self.flowsheet().config.time,
            initialize=0.9,
            doc="Isentropic efficiency"
        )
        eff.fix()

        pratio = self.ratioP = pyo.Var(
            self.flowsheet().config.time,
            initialize=0.7,
            doc="Ratio of outlet to inlet pressure"
        )

        # Some shorter refernces to property blocks
        properties_in = self.control_volume.properties_in
        properties_out = self.control_volume.properties_out

for e in entity_sensors.keys():
            entity_sensors[e] = list(sorted(entity_sensors[e]))

        model.entity_list = pe.Set(initialize=entity_list, ordered=True)
        model.sensor_list = pe.Set(initialize=sensor_list, ordered=True)

        if redundancy > 0:
            model.x = pe.Var(model.entity_list, within=pe.Binary)
        else:
            model.x = pe.Var(model.entity_list, bounds=(0,1))
        model.y = pe.Var(model.sensor_list, within=pe.Binary)

        if use_entity_weight:
            entity_weights = entity.set_index('Entity')['Weight']
            model.obj = pe.Objective(expr=sum(float(entity_weights[e])*model.x[e] for e in entity_list), sense=pe.maximize)
        else:
            model.obj = pe.Objective(expr=sum(model.x[e] for e in entity_list), sense=pe.maximize)

        def entity_covered_rule(m, e):
            if redundancy > 0:
                return (redundancy + 1.0)*m.x[e] <= sum(m.y[b] for b in entity_sensors[e])
            return m.x[e] <= sum(m.y[b] for b in entity_sensors[e])
        model.entity_covered = pe.Constraint(model.entity_list, rule=entity_covered_rule)

        if sensor_budget is None:
            if use_sensor_cost:
                raise ValueError('CoverageFormulation: sensor_budget must be specified if use_sensor_cost is set to True.')
            sensor_budget = len(sensor_list) # no sensor budget provided - allow all sensors
        model.sensor_budget = pe.Param(initialize=sensor_budget, mutable=True)

        if use_sensor_cost:

s_entities = coverage_series[s]

            for e in s_entities:
                entity_sensors[e].add(s)

        for e in entity_sensors.keys():
            entity_sensors[e] = list(sorted(entity_sensors[e]))

        model.entity_list = pe.Set(initialize=entity_list, ordered=True)
        model.sensor_list = pe.Set(initialize=sensor_list, ordered=True)

        if redundancy > 0:
            model.x = pe.Var(model.entity_list, within=pe.Binary)
        else:
            model.x = pe.Var(model.entity_list, bounds=(0,1))
        model.y = pe.Var(model.sensor_list, within=pe.Binary)

        if use_entity_weight:
            entity_weights = entity.set_index('Entity')['Weight']
            model.obj = pe.Objective(expr=sum(float(entity_weights[e])*model.x[e] for e in entity_list), sense=pe.maximize)
        else:
            model.obj = pe.Objective(expr=sum(model.x[e] for e in entity_list), sense=pe.maximize)

        def entity_covered_rule(m, e):
            if redundancy > 0:
                return (redundancy + 1.0)*m.x[e] <= sum(m.y[b] for b in entity_sensors[e])
            return m.x[e] <= sum(m.y[b] for b in entity_sensors[e])
        model.entity_covered = pe.Constraint(model.entity_list, rule=entity_covered_rule)

        if sensor_budget is None:
            if use_sensor_cost:
                raise ValueError('CoverageFormulation: sensor_budget must be specified if use_sensor_cost is set to True.')

def _calculate_scale_factors_from_nominal(m):
    """PRIVATE FUNCTION
    For variables and expressions, if a nominal value is provided calculate the
    scaling factor.

    Args:
        m (Block): a pyomo block to calculate scaling factors for

    Returns:
        None
    """

    for c in m.component_data_objects((pyo.Var, pyo.Expression, pyo.Objective)):
        # Check for a scaling expression.  If there is one, use it to calculate
        # a scaling factor otherwise use autoscaling.
        if not hasattr(c.parent_block(), "nominal_value"):
            continue  # no scaling expression supplied
        elif c not in c.parent_block().nominal_value:
            continue  # no scaling expression supplied

        if not hasattr(c.parent_block(), "scaling_factor"):
            # if there is no scaling_factor Suffix yet make one
            c.parent_block().scaling_factor = pyo.Suffix(
                direction=pyo.Suffix.EXPORT)

        # Add scaling factor from nominal value of variables or expressions
        c.parent_block().scaling_factor[c] = 1 / c.parent_block().nominal_value[
            c]