Co-optimisation of Pyomo and EnergyPlus

Question

I am trying to run pyomo optimisation while considering energyplus(eplus) temperature simulations as constraints of my objective function. It appears that eplus is called only once. How can I ensure that eplus values are updated during iterations?

Below is the code;

import pyomo.environ as pyo
import numpy as np
import matplotlib.pyplot as plt
from EPlus.model import TemperatureModel
from pyomo.opt import SolverFactory

# Initialize the EnergyPlus model
EPLUS_model = TemperatureModel("C:\EnergyPlusV9-4-0")

# Create a Pyomo ConcreteModel
model = pyo.ConcreteModel()

# The time
model.time = pyo.Set(initialize=(range(30, 31 + 1, 30)))

# Optimised Variable
model.power_w = pyo.Var(model.time, within=pyo.NonNegativeReals, bounds=(0, 40000), initialize=5000)

# Defining the expression for indoor temperatures
def ep_model(model, t):
    return EPLUS_model.T(np.zeros(96)+model.power_w[t])[0]
model.indoor_temps = pyo.Expression(model.time, rule=lambda model, t: ep_model(model, t))

# Adding temp limit as a constraint
def LB_rule(model, t):  
    return model.indoor_temps[t] >= 13
model.temp_constraint = pyo.Constraint(model.time, rule=LB_rule)

# Define the objective to minimize total power consumption
def obj_rule(model):
    return sum(model.power_w[t] for t in model.time)

model.total_power = pyo.Objective(rule=obj_rule, sense=pyo.minimize)

solver = SolverFactory('ipopt')
result = solver.solve(model, tee=False)

solved_p_values = {t: model.power_w[t].value for t in model.time}

solved_t_values = solved_t_values = [pyo.value(model.indoor_temps[t]) for t in model.time]#{t: model.indoor_temps[t]() for t in model.time}