scipy.optimize functions take an args parameter, as a way of passing extra arguments to the function. You could use that to pass a number of the parameters to heat_balance.
For example, let's assume r, T_a and dx are more 'variable' than R_w; that is, more likely to vary from run to run or test case. I'll also switch the order of T_1 and T_2, since the later is, apparently, the variable you want to minimize.
def heat_balance(T_2,T_1, T_a, r, dx):
heat_in = (T_1-T_2) * scipy.pi * r**2 * R_w
heat_out = ( (T_1 + T_2) / 2 - T_a ) * 2 * r * scipy.pi * dx * R_f
return abs(heat_in-heat_out)
Then in simulation minimize could be called with:
result = optimize.minimize(heat_balance, T[-1]-dx, (T[-1], T_a, r, dx),
method="SLSQP")
newT = result['x'] # not result.x?
# what's the purpose of the float() call?
T.append(newT)
The simulation definition could also include these 'global' variables.
If we are going to define object classes, I'd think more in terms of 'pipe' objects rather than 'simulation' objects. A 'pipe' would have certain properties - dimensions, conductivity etc, that could differ from one pipe to another. So it could be convenient to define several such objects, and pass them to the simulation.
I'd rewrite the simulation iteration as something like:
def simulation(...):
X = np.arange(dx, L, dx)
T = np.zeros_like(X)
for i, x in enumerate(X):
result = optimize.minimize(heat_balance, T[-1]-dx, (T[-1], T_a, r, dx), method="SLSQP")
newT = result['x']
T[i] = newT
return T,X # return a simple tuple
And later call:
T,X = simulation(<parameters>)
plot(T, X,...)
I haven't tested these changes, but I think they give an adequate idea of how the code and calling structure could be cleaned up.
These changes don't address the speed issue. Repeatedly calling minimize over the length of pipe, does sound expensive.
Here's a version of your script that could make it easier to play with some of the parameters, such as the dx. I can import it into a shell (e.g. ipython session) and perform timings on pieces. I'm experimenting with 2 forms of the func. Precalculating values may say 20-30% time, not a change. dx has most effect on simulation times.
I also collect the nfev statistic, incase that gives any ideas of how costly the minimize call is. On average it calls the balance function 25 times per call, or 6000 times for a dx=.1.
Fine tuning the heat_balance equation helps some with speed, but I think being smarter with dx and the use of minimize help more.
import numpy as np
from scipy.optimize import minimize
from matplotlib import pyplot as plt
Rw = 0.654 # W/mk water heat loss ocefficient
Rf = 0.042 # W/mk pipe heat loss coefficient
r = 0.2 # m pipe radius
T0 = 50 # C initial temperature
Ta = 20 # C ambiant temperature
dx = 0.01 # m differential distance
L = 25 # m length of pipe
pi = np.pi
dt = 0 # tweak to last T
def heat_balance(T2, T1, Rf, r, Ta):
heat_in = (T1-T2) * pi * r**2 * Rw
heat_out = ( (T1 + T2) / 2 - Ta ) * 2 * r * pi * dx * Rf
return abs(heat_in-heat_out)
class Balance(object):
# heat calc with precalc
# seems to save 30% on time
def __init__(self, Rf, r, Ta):
self.coeff1 = pi * r**2 * Rw
self.coeff2 = 2 * r * pi * dx * Rf
self.Ta = Ta
def __call__(self, s, T1):
heat_in = (T1-s) * self.coeff1
heat_out = ( (T1+s)/2 - self.Ta) * self.coeff2
return abs(heat_in - heat_out)
def simulation1(Rf, T0, r, Ta):
# simulation with plain function
X = np.arange(0, L, dx)
T = np.zeros_like(X)
T[0] = lastT = T0
cnt = np.zeros_like(X)
func = heat_balance
for i,x in enumerate(X[1:]):
args = (lastT, Rf, r, Ta)
result = minimize(func, lastT-dt, args, method="SLSQP")
T[i+1] = lastT = result.x
cnt[i+1] = result.nfev
cnt = cnt[1:]
return X, T, cnt
def simulation2(Rf, T0, r, Ta):
# simulation with Balance class
X = np.arange(0, L, dx)
T = np.zeros_like(X)
T[0] = lastT = T0
cnt = np.zeros_like(X)
func = Balance(Rf, r, Ta)
for i,x in enumerate(X[1:]):
args = (lastT,)
result = minimize(func, lastT-dt, args, method="SLSQP")
T[i+1] = lastT = result.x
cnt[i+1] = result.nfev
cnt = cnt[1:]
return X, T, cnt
simulation = simulation2
if __name__ == '__main__':
import sys
if sys.argv[1:]:
dx = float(sys.argv[1])
if sys.argv[2:]:
dt = float(sys.argv[2]) # temp offset
Rfs = [7,14,50,100,200]
typs = ['+','x','|','_','*']
cols = ['red','blue','black','purple','brown']
for rf, typ, col in zip(Rfs, typs, cols):
Rf = 1.0/rf
print(Rf)
X, T, cnt = simulation(Rf, T0, r, Ta)
print('nfev count: ',len(cnt), cnt.min(), cnt.max(), cnt.sum(), cnt.mean())
plt.scatter(X, T, c=col, marker=typ, s = 100, label=rf, lw=3)
plt.legend(fontsize=20,title='Pipe Insolation R Value')
plt.xlabel('Distance (m)',fontsize=20)
plt.ylabel('Average Temperature (C)',fontsize=20)
plt.show()
