1import sys
2import pyomo.environ as pyo
3import numpy.random as rnd
4import pandas as pd
5import pyomo.contrib.pynumero.examples.external_grey_box.param_est.models as po
6
7def perform_estimation_pyomo_only(data_fname, solver_trace=False):
8    # read in our data file - careful with formats
9    df = pd.read_csv(data_fname)
10    npts = len(df)
11
12    # create our parameter estimation formulation
13    m = pyo.ConcreteModel()
14    m.df = df
15    m.PTS = pyo.Set(initialize=range(npts), ordered=True)
16
17    # create a separate Pyomo block for each data point
18    def _model_i(b, i):
19        po.build_single_point_model_pyomo_only(b)
20    m.model_i = pyo.Block(m.PTS, rule=_model_i)
21
22    # we want the parameters to be the same across all the data pts
23    m.UA = pyo.Var()
24    def _eq_parameter(m, i):
25        return m.UA == m.model_i[i].UA
26    m.eq_parameter = pyo.Constraint(m.PTS, rule=_eq_parameter)
27
28    # define the least squares objective function
29    def _least_squares(m):
30        obj = 0
31        for i in m.PTS:
32            row = m.df.iloc[i]
33
34            # error in inputs measured
35            obj += (m.model_i[i].Th_in - float(row['Th_in']))**2
36            obj += (m.model_i[i].Tc_in - float(row['Tc_in']))**2
37
38            # error in outputs
39            obj += (m.model_i[i].Th_out - float(row['Th_out']))**2
40            obj += (m.model_i[i].Tc_out - float(row['Tc_out']))**2
41        return obj
42    m.obj = pyo.Objective(rule=_least_squares)
43
44    solver = pyo.SolverFactory('ipopt')
45    status = solver.solve(m, tee=solver_trace)
46
47    return m
48
49def perform_estimation_external(data_fname, solver_trace=False):
50    # read in our data file - careful with formats
51    df = pd.read_csv(data_fname)
52    npts = len(df)
53
54    # create our parameter estimation formulation
55    m = pyo.ConcreteModel()
56    m.df = df
57    m.PTS = pyo.Set(initialize=range(npts), ordered=True)
58
59    # create a separate Pyomo block for each data point
60    def _model_i(b, i):
61        po.build_single_point_model_external(b)
62    m.model_i = pyo.Block(m.PTS, rule=_model_i)
63
64    # we want the parameters to be the same across all the data pts
65    # create a global parameter and provide equality constraints to
66    # the parameters in each model instance
67    m.UA = pyo.Var()
68    def _eq_parameter(m, i):
69        return m.UA == m.model_i[i].egb.inputs['UA']
70    m.eq_parameter = pyo.Constraint(m.PTS, rule=_eq_parameter)
71
72    # define the least squares objective function
73    def _least_squares(m):
74        obj = 0
75        for i in m.PTS:
76            row = m.df.iloc[i]
77
78            # error in inputs measured
79            obj += (m.model_i[i].egb.inputs['Th_in'] - float(row['Th_in']))**2
80            obj += (m.model_i[i].egb.inputs['Tc_in'] - float(row['Tc_in']))**2
81
82            # error in outputs
83            obj += (m.model_i[i].egb.inputs['Th_out'] - float(row['Th_out']))**2
84            obj += (m.model_i[i].egb.inputs['Tc_out'] - float(row['Tc_out']))**2
85        return obj
86    m.obj = pyo.Objective(rule=_least_squares)
87
88    solver = pyo.SolverFactory('cyipopt')
89    status, nlp = solver.solve(m, tee=solver_trace, return_nlp=True)
90
91    if solver_trace:
92        # use the NLP object to access additional information if so desired
93        # for example:
94        names = nlp.primals_names()
95        values = nlp.evaluate_grad_objective()
96        print({names[i]:values[i] for i in range(len(names))})
97
98    return m
99
100if __name__ == '__main__':
101    m = perform_estimation_pyomo_only(sys.argv[1])
102    print(pyo.value(m.UA))
103    m = perform_estimation_external(sys.argv[1])
104    print(pyo.value(m.UA))
105
106