ODE's not being solved because of a data mismatch in python. I have tried several approaches but keep getting the same error, what's wrong?

I am trying to use python to solve a set of ode’s for an electrical engineering project I am working on. In this project I am trying to solve a set of ODE’s in order to get the transfomer hotspot temperature against time for a day. I have a data frame showing the hourly total power consumed by a particular load.

import numpy as np
import pandas as pd
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Constants from the problem
transformer_size = 500
x = 0.8
y = 1.6
R = 8
k11 = 1
k21 = 1
k22 = 2
tau_w = 4
tau_o = 180
delta_theta_or = 45
delta_theta_hr = 45
theta_alpha = 30  # Ambient temperature (°C)

# Assuming 'total_load_df' is already defined
# Calculating the load factor K
K = total_load_df['Total_PT'] / transformer_size

# Prepare DataFrame for K values
K_df = pd.DataFrame({
    'hour': range(24),
    'K': K  # The load factor K over 24 hours
})

# Debugging: Check if K_df is populated correctly
print("K_df values:\n", K_df)

def thermal_model(t, y, K_df):
    theta_o, Delta_theta_h1, Delta_theta_h2 = y

    # Interpolating K for the current time t
    K_t = np.interp(t, K_df['hour'], K_df['K'])
    
    # Ensure K_t is a scalar
    K_t = K_t.item() if isinstance(K_t, np.ndarray) else K_t
    print(f"At time t={t}, interpolated K_t={K_t}")  # Debugging: Check interpolated K_t

    # Top-oil temperature differential equation (scalar)
    dtheta_o_dt = (1 + K_t**2 * R / (1 + R))**x * delta_theta_or / (k11 * tau_o) - (theta_o - theta_alpha) / tau_o

    # Hot-spot temperature rise equations (scalars)
    dDelta_theta_h1_dt = (k21 * K_t**y * delta_theta_hr - Delta_theta_h1) / (k22 * tau_w)
    dDelta_theta_h2_dt = ((k21 - 1) * K_t**y * delta_theta_hr - Delta_theta_h2) / (tau_o * k22)

    # Debugging: Ensure the outputs are scalar values
    print(f"At time t={t}, dtheta_o_dt={dtheta_o_dt}, dDelta_theta_h1_dt={dDelta_theta_h1_dt}, dDelta_theta_h2_dt={dDelta_theta_h2_dt}")
    
    # Return as a list of scalars
    return [dtheta_o_dt, dDelta_theta_h1_dt, dDelta_theta_h2_dt]


# Initial conditions
theta_o_init = 60  # Initial top-oil temperature (°C)
Delta_theta_h1_init = 20  # Initial hot-spot temperature rise (part 1)
Delta_theta_h2_init = 10  # Initial hot-spot temperature rise (part 2)

# Time range for 24 hours
t_span = [0, 24]
t_eval = np.linspace(0, 24, 1000)  # Time points where the solution is evaluated

# Solve the ODEs
sol = solve_ivp(
    thermal_model, t_span, [theta_o_init, Delta_theta_h1_init, Delta_theta_h2_init], 
    t_eval=t_eval, args=(K_df,)
)


# Extract the solution
theta_o_sol = sol.y[0]  # Top-oil temperature
Delta_theta_h1_sol = sol.y[1]  # Hot-spot temperature rise part 1
Delta_theta_h2_sol = sol.y[2]  # Hot-spot temperature rise part 2

# Calculate total hot-spot temperature and hotspot temperature rise
Delta_theta_h_sol = Delta_theta_h1_sol - Delta_theta_h2_sol
theta_h_sol = theta_o_sol + Delta_theta_h_sol

# Plot results
plt.figure(figsize=(10, 6))
plt.plot(t_eval, theta_o_sol, label="Top-oil Temperature (θ_o,t)", color="blue")
plt.plot(t_eval, theta_h_sol, label="Hot-spot Temperature (θ_h,t)", color="red")
plt.xlabel('Time (hours)')
plt.ylabel('Temperature (°C)')
plt.title('Top-oil and Hot-spot Temperature over 24 hours')
plt.legend()
plt.grid(True)
plt.show()

When I ran my code I am getting the following error message:

K_df values:
     hour         K
0      0  0.248627
1      1  0.195822
2      2  0.155273
3      3  0.140915
4      4  0.134600
5      5  0.169652
6      6  0.346817
7      7  0.645293
8      8  0.742582
9      9  0.543450
10    10  0.451147
11    11  0.382673
12    12  0.405557
13    13  0.486433
14    14  0.518242
15    15  0.357977
16    16  0.369747
17    17  0.438374
18    18  0.591805
19    19  0.692828
20    20  0.619413
21    21  0.628065
22    22  0.447500
23    23  0.331827
At time t=0.0, interpolated K_t=0.24862719825155716
At time t=0.0, dtheta_o_dt=0.09426365740423229, dDelta_theta_h1_dt=[-2.5        -2.5        -2.49999492], dDelta_theta_h2_dt=[-0.02777778 -0.02777778 -0.02777778]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-22-ef76ccdba745> in <cell line: 67>()
     65 
     66 # Solve the ODEs
---> 67 sol = solve_ivp(
     68     thermal_model, t_span, [theta_o_init, Delta_theta_h1_init, Delta_theta_h2_init],
     69     t_eval=t_eval, args=(K_df,)

3 frames
/usr/local/lib/python3.10/dist-packages/scipy/integrate/_ivp/base.py in fun_wrapped(t, y)
     21 
     22     def fun_wrapped(t, y):
---> 23         return np.asarray(fun(t, y), dtype=dtype)
     24 
     25     return fun_wrapped, y0

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (3,) + inhomogeneous part.

The function

thermal_model

is returning an array with inconsistent shapes, specifically after the first dimension.

so this was tried, thought it didn’t solve the error

def thermal_model(t, y, K_df):
    theta_o, Delta_theta_h1, Delta_theta_h2 = y

    # Interpolating K for the current time t
    K_t = np.interp(t, K_df['hour'], K_df['K'])
    
    # Ensure K_t is a scalar by taking the first element if it's an array:
    K_t = K_t.item() if isinstance(K_t, np.ndarray) else K_t
    
    # Top-oil temperature differential equation (scalar)
    dtheta_o_dt = (1 + K_t**2 * R / (1 + R))**x * delta_theta_or / (k11 * tau_o) - (theta_o - theta_alpha) / tau_o

    # Hot-spot temperature rise equations (scalars)
    dDelta_theta_h1_dt = (k21 * K_t**y * delta_theta_hr - Delta_theta_h1) / (k22 * tau_w)
    dDelta_theta_h2_dt = ((k21 - 1) * K_t**y * delta_theta_hr - Delta_theta_h2) / (tau_o * k22)

    # Return as a list of scalars
    return [dtheta_o_dt, dDelta_theta_h1_dt, dDelta_theta_h2_dt]

You need to sign in to view this answers

Related Post