If2 / 3 using the Orthogonal Collocation on Finite Elements method in Gekko

I am currently working on a dynamic heat transfer model, using the Gekko library. The model consists of a system of partial differential equations. As indicated on the APMonitor platform, the differentials that vary in space have been discretized, while differentials in time have been solved numerically by Orthogonal Collocation on Finite Elements method.

I have two problems:

1. The model solves only with a maximum of 3 discretized points in space and 30 in time with NODES = 3. While for NODES = 4 the nodal points do not follow the values of the function. What explains this type of behavior?

2. The model is not solved when any of the if2/3 are used. All types of solvers were tested and none found the solution. When m.if3 are replaced by m.Intermediates the model solves normally. Furthermore, if the number of nodes is not specified, the model is solved with up to n = 20 and nt = 300 points. What can explain this drop in performance by adding ifs and nodal points? Is the model incorrectly formulated or structured?

What can be done to improve? The model needs to be solved with at least 100 points in time and NODES = 3.

***I would also like to know how to use an if2/3 inside another if2/3.

This is the main part of the model:

import numpy as np
import data
from func_H import func_h
from func_KW import func_kw
from func_den import func_fg_cpv
from func_PG import func_pg
from func_CPV import func_cpv11, func_cpv21, func_cpv31, func_cpv41, func_cpv51, func_cpv6, func_cpv7, func_cpv8, func_cpv12, func_cpv22, func_cpv32, func_cpv42, func_cpv52
from func_CPL import func_cpl11, func_cpl12, func_cpl21, func_cpl22, func_cpl31, func_cpl32, func_cpl4, func_cpl5, func_cpl6, func_cpl71, func_cpl72, func_cpl8, func_cpl9
from func_cond_in import func_cond_inic
from func_DHW import func_dhw

import matplotlib.pyplot as plt

from gekko import GEKKO
m = GEKKO(remote=False)

zf = 0.17
tf = 768.24
dz = zf/n

n = 3
nt = 10
m.time = np.linspace(0,tf,nt)
nd = 3
####################################################################################################

u = func_cond_inic(data.tg,data.g,data.p1,data.yi[0],data.yi[1],data.yi[2],data.yi[3],data.yi[4],data.yi[5],data.yi[6],data.yi[7],n)

tg = [m.Var(u[0][i]) for i in range(n)]
g = [m.Var(u[1][i]) for i in range(n)]
pg = [m.Var(u[2][i]) for i in range(n)]

tp = [m.Var(data.tp) for i in range(n)]

y1 = [m.Var(u[3][i][0]) for i in range(n)]
y2 = [m.Var(u[4][i][0]) for i in range(n)]
y3 = [m.Var(u[5][i][0]) for i in range(n)]
y4 = [m.Var(u[6][i][0]) for i in range(n)]
y5 = [m.Var(u[7][i][0]) for i in range(n)]
y6 = [m.Var(u[8][i][0]) for i in range(n)]
y7 = [m.Var(u[9][i][0]) for i in range(n)]
y8 = [m.Var(u[10][i][0]) for i in range(n)]

y = [np.array([y1[i], y2[i], y3[i], y4[i], y5[i], y6[i], y7[i], y8[i]]) for i in range(n)]
yt = [np.transpose(y[i]) for i in range(n)]
    
x1 = [m.Var(data.xi[0]) for i in range(n)]
x2 = [m.Var(data.xi[1]) for i in range(n)]
x3 = [m.Var(data.xi[2]) for i in range(n)]
x4 = [m.Var(data.xi[3]) for i in range(n)]
x5 = [m.Var(data.xi[4]) for i in range(n)]
x6 = [m.Var(data.xi[5]) for i in range(n)]
x7 = [m.Var(data.xi[6]) for i in range(n)]
x8 = [m.Var(data.xi[7]) for i in range(n)]
x9 = [m.Var(data.xi[8]) for i in range(n)]

x = [np.array([x1[i], x2[i], x3[i], x4[i], x5[i], x6[i], x7[i], x8[i], x9[i]]) for i in range(n)]

####################################################################################################

rhog = [m.Intermediate((pg[i]*(data.mmg @ y[i])/(data.r*tg[i]))) for i in range(n)]

kw = [m.Intermediate(func_kw(g[i],tg[i],y[i],pg[i],rhog[i])) for i in range(n)]

weg = [m.Intermediate(m.exp(-3.5+(0.0549/tg[i]))) for i in range(n)]

wg = [m.Intermediate((pg[i]*data.mmg[3]*y4[i])/(data.r*tg[i])) for i in range(n)]

rw_f = [m.Intermediate(kw[i]*data.a*(weg[i]-wg[i])) for i in range(n)]

rw_m = [m.Intermediate(rw_f[i]*1000/18.01528) for i in range(n)]

dhw = [m.Intermediate(func_dhw(tg[i],tp[i])) for i in range(n)]

qcond = [m.Intermediate(rw_m[i]*dhw[i]) for i in range(n)]

# # Heat capacity 
cpv1 = [m.if2(tg[i]-700,func_cpv11(tg[i]),func_cpv12(tg[i])) for i in range(n)]
cpv2 = [m.if2(tg[i]-500,func_cpv21(tg[i]),func_cpv22(tg[i])) for i in range(n)]
cpv3 = [m.if2(tg[i]-1200,func_cpv31(tg[i]),func_cpv32(tg[i])) for i in range(n)]
cpv4 = [m.if2(tg[i]-1700,func_cpv41(tg[i]),func_cpv42(tg[i])) for i in range(n)]
cpv5 = [m.if2(tg[i]-1300,func_cpv51(tg[i]),func_cpv52(tg[i])) for i in range(n)]
# cpv1 = [m.Intermediate(func_cpv11(tg[i])) for i in range(n)]    
# cpv2 = [m.Intermediate(func_cpv21(tg[i])) for i in range(n)]
# cpv3 = [m.Intermediate(func_cpv31(tg[i])) for i in range(n)]
# cpv4 = [m.Intermediate(func_cpv41(tg[i])) for i in range(n)]
# cpv5 = [m.Intermediate(func_cpv51(tg[i])) for i in range(n)]
cpv6 = [m.Intermediate(func_cpv6(tg[i])) for i in range(n)]
cpv7 = [m.Intermediate(func_cpv7(tg[i])) for i in range(n)]
cpv8 = [m.Intermediate(func_cpv8(tg[i])) for i in range(n)]

cpv = [[cpv1[i], cpv2[i], cpv3[i], cpv4[i], cpv5[i], cpv6[i], cpv7[i], cpv8[i]] for i in range(n)]

fg_cpv = [m.Intermediate(func_fg_cpv(g[i],cpv[i],y[i])) for i in range(n)]

h = [m.Intermediate(func_h(g[i],tg[i],y[i],cpv[i])) for i in range(n)]

qconv = [m.Intermediate((data.a*h[i]*(tg[i]-tp[i]))) for i in range(n)]

p = [m.Intermediate(func_pg(g[i],tg[i],y[i],rhog[i])) for i in range(n)]

# # Heat capacity 
cpl1 = [m.if2(tp[i]-950,func_cpl11(tp[i]),func_cpl12(tp[i])) for i in range(n)]    
cpl2 = [m.if2(tp[i]-900,func_cpl21(tp[i]),func_cpl22(tp[i])) for i in range(n)]
cpl3 = [m.if2(tp[i]-847,func_cpl31(tp[i]),func_cpl32(tp[i])) for i in range(n)]
# cpl1 = [m.Intermediate(func_cpl11(tp[i])) for i in range(n)]
# cpl2 = [m.Intermediate(func_cpl21(tp[i])) for i in range(n)]
# cpl3 = [m.Intermediate(func_cpl31(tp[i])) for i in range(n)]
cpl4 = [m.Intermediate(func_cpl4(tp[i])) for i in range(n)]
cpl5 = [m.Intermediate(func_cpl5(tp[i])) for i in range(n)]
cpl6 = [m.Intermediate(func_cpl6(tp[i])) for i in range(n)]
cpl7 = [m.if2(tp[i]-500,func_cpl71(tp[i]),func_cpl72(tp[i])) for i in range(n)]
# cpl7 = [m.Intermediate(func_cpl71(tp[i])) for i in range(n)]
cpl8 = [m.Intermediate(func_cpl8(tp[i])) for i in range(n)]
cpl9 = [m.Intermediate(func_cpl9(tp[i])) for i in range(n)]

cpl = [[cpl1[i], cpl2[i], cpl3[i], cpl4[i], cpl5[i], cpl6[i], cpl7[i], cpl8[i], cpl9[i]] for i in range(n)]

rhop1 = [m.Intermediate((x1[i]*data.rhoc[0]*(1-data.e))) for i in range(n)]
rhop2 = [m.Intermediate((x2[i]*data.rhoc[1]*(1-data.e))) for i in range(n)]
rhop3 = [m.Intermediate((x3[i]*data.rhoc[2]*(1-data.e))) for i in range(n)]
rhop4 = [m.Intermediate((x4[i]*data.rhoc[3]*(1-data.e))) for i in range(n)]
rhop5 = [m.Intermediate((x5[i]*data.rhoc[4]*(1-data.e))) for i in range(n)]
rhop6 = [m.Intermediate((x6[i]*data.rhoc[5]*(1-data.e))) for i in range(n)]
rhop7 = [m.Intermediate((x7[i]*data.rhoc[6]*(1-data.e))) for i in range(n)]
rhop8 = [m.Intermediate((x8[i]*data.rhoc[7]*(1-data.e))) for i in range(n)]
rhop9 = [m.Intermediate((x9[i]*data.rhoc[8]*(1-data.e))) for i in range(n)]
rhop = [np.array([rhop1[i], rhop2[i], rhop3[i], rhop4[i], rhop5[i], rhop6[i], rhop7[i], rhop8[i], rhop9[i]]) for i in range(n)]

rho_cpl = [m.Intermediate((rhop[i] @ cpl[i])) for i in range(n)]

####################################################################################################

m.Equation(tg[0] == data.tg)
m.Equation([tg[i] == tg[i-1] + dz*(qcond[i-1]-qconv[i-1])/fg_cpv[i-1] for i in range(1,n)])

m.Equation(g[0] == data.g)
m.Equation([g[i] == g[i-1] + dz*rw_f[i-1] for i in range(1,n)])

m.Equation(pg[0] == data.p1)
m.Equation([pg[i] == pg[i-1] + dz*p[i-1] for i in range(1,n)])
    
m.Equation(y1[0] == data.yi[0])
m.Equation(y2[0] == data.yi[1])
m.Equation(y3[0] == data.yi[2])
m.Equation(y4[0] == data.yi[3])
m.Equation(y5[0] == data.yi[4])
m.Equation(y6[0] == data.yi[5])
m.Equation(y7[0] == data.yi[6])
m.Equation(y8[0] == data.yi[7])
m.Equation([y1[i] == y1[i-1] - dz*rw_f[i-1]*y1[i-1]/g[i-1] for i in range(1,n)])
m.Equation([y2[i] == y2[i-1] - dz*rw_f[i-1]*y2[i-1]/g[i-1] for i in range(1,n)])
m.Equation([y3[i] == y3[i-1] - dz*rw_f[i-1]*y3[i-1]/g[i-1] for i in range(1,n)])
m.Equation([y4[i] == y4[i-1] + dz*rw_f[i-1]*(1-y4[i-1])/g[i-1] for i in range(1,n)])
m.Equation([y5[i] == y5[i-1] - dz*rw_f[i-1]*y5[i-1]/g[i-1] for i in range(1,n)])
m.Equation([y6[i] == y6[i-1] - dz*rw_f[i-1]*y6[i-1]/g[i-1] for i in range(1,n)])
m.Equation([y7[i] == y7[i-1] - dz*rw_f[i-1]*y7[i-1]/g[i-1] for i in range(1,n)])
m.Equation([y8[i] == y8[i-1] - dz*rw_f[i-1]*y8[i-1]/g[i-1] for i in range(1,n)])

m.Equation([tp[i].dt() == ((qconv[i]-qcond[i])/rho_cpl[i]) for i in range(n)])

m.Equation([x1[i].dt() == 0 for i in range(n)])
m.Equation([x2[i].dt() == 0 for i in range(n)])
m.Equation([x3[i].dt() == 0 for i in range(n)])
m.Equation([x4[i].dt() == 0 for i in range(n)])
m.Equation([x5[i].dt() == 0 for i in range(n)])
m.Equation([x6[i].dt() == 0 for i in range(n)])
m.Equation([x7[i].dt() == (-rw_f[i]/(rhop1[i]+rhop2[i]+rhop3[i]+rhop4[i]+rhop5[i]+rhop6[i]+rhop7[i]+rhop8[i]+rhop9[i])) for i in range(n)])
m.Equation([x8[i].dt() == 0 for i in range(n)])
m.Equation([x9[i].dt() == 0 for i in range(n)])

####################################################################################################
m.options.imode = 4
m.options.SOLVER = 1
m.options.NODES = nd
m.options.REDUCE = 2
m.options.MAX_ITER = 250
m.solve()

The complete model is available at the link:

https://drive.google.com/drive/folders/1_AFnBPZ_mcEhwZb80a3VYOcg8PVXP5Qj?usp=sharing

Thank you in advance.