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