PYOMO ODE Solve with Manual Discretization

[Audrey Demafo]

Hi, can anyone help me figure out where I made a mistake in the following code ? I am trying to solve a simple ode (dxdt = 5x -3, with initial conditions x(0)=0.2). However, the solution I am getting is wrong as it does not match the values obtained with the analytical solution which (x(t) = -0.4*exp(5*t) + 0.6). 

m = ConcreteModel()
m.tf = Param(initialize=24)
m.t = RangeSet(0, value(m.tf))
m.x = Var(m.t)
m.ht = Param(initialize=(value(m.tf))/(value(m.tf)))

m.dxdt = Var(m.t)

## Define the ODE
def _ode(m, k):
    return m.dxdt[k]==5*m.x[k] - 3

## Discretization scheme
def _ode_discr(m, k):
    if k==0:
        return Constraint.Skip
    return m.x[k] == m.x[k-1] + m.ht*m.dxdt[k-1]   ## Forward euler discretization

def _initial_cond(m):
    return m.x[0]==0.2  
m.initial_cond = Constraint(rule=_initial_cond)   

m.obj = 1

sol=[]

for t in m.t:
    if 'ode' and 'ode_discr' in m.component_map():
        m.del_component(m.ode)
        m.del_component(m.ode_discr)
    m.ode = Constraint(m.t, rule=_ode)
    m.ode_discr = Constraint(m.t, rule=_ode_discr)
    solver = SolverFactory("ipopt")
    solver.solve(m, tee=True)
    print("at", t, "", "value=", value(m.x[t]))
    sol.append(value(m.x[t]))
print(sol)

Any help on this would be greatly appreciated. 

Best regards,

Audrey.

DISCLAIMER: The contents of this email and any attachments are confidential. They are intended for the named recipient(s) only. If you have received this email by mistake, please notify the sender immediately and you are herewith notified that the contents are legally privileged and that you do not have permission to disclose the contents to anyone, make copies thereof, retain or distribute or act upon it by any means, electronically, digitally or in print. The views expressed in this communication may be of a personal nature and not be representative of AIMS-Rwanda and/or any of its Initiatives.

[Ahmad Heidari]

Hello Audrey,

I hope you are doing well.

Regarding the nature of your problem, you could decrease the value of "ht" for better stability as follows. Please place all the code in one cell. 

from pyomo.environ import *
import matplotlib.pyplot as plt

m = ConcreteModel()

#Parameters

m.tf = Param(initialize=24)
m.ht = Param(initialize=0.01)  # smaller step for better stability
m.t = RangeSet(0, int(value(m.tf)/value(m.ht)))

#Variables

m.x    = Var(m.t)
m.dxdt = Var(m.t)

#ODE Discretization (Euler forward)

def _ode_discr(m, k):
    if k == 0:

        return Constraint.Skip
    return m.x[k] == m.x[k-1] + m.ht * m.dxdt[k-1]

m.ode_discr = Constraint(m.t, rule=_ode_discr)

#ODE equation: dx/dt = 5*x - 3

def _ode(m, k):
    return m.dxdt[k] == 5 * m.x[k] - 3

m.ode = Constraint(m.t, rule=_ode)

#Initial condition

m.initial_cond = Constraint(expr=m.x[0] == 0.2)

#Objective (dummy)

m.obj = Objective(expr=sum(m.dxdt[t]**2 for t in m.t), sense=minimize)

#Solve the model

solver = SolverFactory('ipopt')
solver.solve(m, tee=True)

#Extract and print results

sol = [value(m.x[t]) for t in m.t]
for t in m.t:
    print(f"At time {t}: x = {value(m.x[t])}")

#Plot the numerical solution

plt.plot(list(m.t), sol, marker='o', label='Numerical solution')
plt.xlabel("Time (t)")
plt.ylabel("x(t)")
plt.title("ODE Solution Using Pyomo")
plt.grid(True)
plt.legend()
plt.show()