CODEWITHSUNDEEP
pyomo
Mharyoowa
Mharyoowa

Reputation: 3

How to do Interpolation in Pyomo.DAE

I have some time steps values that I already defined in my data but I don't wanna make use of it. I instead want to make use of the discretized points from the ContinuousSet as my new time span.

I tried interpolating my data after applying the discretization but after running it, I noticed three things:

  1. my constraints were more than my variables.
  2. when I print my parameters, it shows the interpolated time points but no value is attached to those time points.
  3. I also noticed the discretization is adding additional points than what I specified (I don't think this is right). I have 97 rows in my data but when I print my parameters, it shows 98.

I can't really figure out where the error is coming from.

Below is my code and a picture of the error:

df = pd.read_csv('digester_influent_PyADM13.csv')
model = m = ConcreteModel()

m.t = ContinuousSet(bounds=(0,1))

m.Ssu_in = Param(m.t, mutable=True)
m.Saa_in = Param(m.t, mutable=True)
m.Sfa_in = Param(m.t, mutable=True)
m.Q = Param(m.t, mutable=True)
m.V_liq = Param(initialize=3400, within=PositiveReals)

discretizer = TransformationFactory('dae.finite_difference')
discretizer.apply_to(m,nfe=(96*1+1),wrt=m.t,scheme='BACKWARD')

# Interpolate the data
timepoints = list(ContinuousSet(bounds=(0,1)))
data_timepoints = df['time'].tolist()
data_profiles1 = df['S_su'].tolist()
data_profiles2 = df['S_aa'].tolist()
data_profiles3 = df['S_fa'].tolist()
data_profiles4 = df['Q'].tolist()

interp_Ssu_values = np.interp(timepoints, data_timepoints, data_profiles1)
interp_Saa_values = np.interp(timepoints, data_timepoints, data_profiles2)
interp_Sfa_values = np.interp(timepoints, data_timepoints, data_profiles3)
interp_Q_values = np.interp(timepoints, data_timepoints, data_profiles4)

for i,t in enumerate(timepoints):
    m.Ssu_in[t] = interp_Ssu_values[i]
    m.Saa_in[t] = interp_Saa_values[i]
    m.Sfa_in[t] = interp_Sfa_values[i]
    m.Q[t] = interp_Q_values[i]

# Variables
m.S_su = Var(m.t, initialize=0.012394, domain=PositiveReals, bounds=(0.001,1))               
m.S_aa = Var(m.t, initialize=0.0055432, domain=PositiveReals, bounds=(0,0.1))                  
m.S_fa = Var(m.t, initialize=0.10741, domain=PositiveReals, bounds=(0.001,2)) 
m.dS_su_dt = DerivativeVar(m.S_su, wrt=m.t)
m.dS_aa_dt = DerivativeVar(m.S_aa, wrt=m.t)
m.dS_fa_dt = DerivativeVar(m.S_fa, wrt=m.t)

# Initial Values
m.S_su[0].fix(0.012394)
m.S_aa[0].fix(0.0055432)
m.S_fa[0].fix(0.10741)

# Differential equations
def S_su_out_bal(m,t):
    return m.dS_su_dt[t] == (m.Q[t]/m.V_liq) * (m.Ssu_in[t] - m.S_su[t]) + 0.000662979
m.Ssu_outcon = Constraint(m.t, rule=S_su_out_bal)

def S_aa_out_bal(m,t):
    return m.dS_aa_dt[t] == (m.Q[t]/m.V_liq) * (m.Saa_in[t] - m.S_aa[t]) - 0.00202160
m.Saa_outcon = Constraint(m.t, rule=S_aa_out_bal)

def S_fa_out_bal(m,t):
    return m.dS_fa_dt[t] == (m.Q[t]/m.V_liq) * (m.Sfa_in[t] - m.S_fa[t]) + 0.005667982
m.Sfa_outcon = Constraint(m.t, rule=S_fa_out_bal)

# Check to see if the number of Constraints equals the number of variables
from pyomo.util.model_size import build_model_size_report
report = build_model_size_report(m)
print('Num constraints: ', report.activated.constraints)
print('Num variables: ', report.activated.variables)

Output:

Num constraints:  294
Num variables:  585

Display Q values:

m.Q.display()

Ouput:

Q : Size=98, Index=t, Domain=Any, Default=None, Mutable=True
    Key      : Value
           0 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.010309 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.020619 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.030928 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.041237 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.051546 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.061856 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.072165 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.082474 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.092784 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.103093 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.113402 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.123711 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.134021 : <class 'pyomo.core.base.param.Param.NoValue'>
     0.14433 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.154639 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.164948 : <class 'pyomo.core.base.param.Param.NoValue'>
    0.175258 : <class 'pyomo.core.base.param.Param.NoValue'>

Upvotes: 0

Views: 106

Answers (1)

Bethany Nicholson
Bethany Nicholson

Reputation: 2828

Here are the issues I see:

  1. This line: timepoints = list(ContinuousSet(bounds=(0,1))) is creating a new ContinuousSet, I think you want this instead: timepoints = list(m.t)

  2. You need to declare all of the Var and DerivativeVar components before applying the discretization transformation.

Upvotes: 1

Related Questions