Hi Jack!
I'm glad to read your comment.
About considering your system as SISO, I don't know your exactly your setup, thus I don't know how weak is that assumption. But I understand that when you say that you are going to hold your power level, you're controlling the hydrogen (or whatever) flow rate so the reaction gives you always the same power, and so, the only temperature variation comes from your air pump.
If you just hold the fuel rate constant, the reacting fuel will increase as you increase the amount of oxygen in the system (probably you already know this, but it is better if I don't assume).
Now, about your question itself, there are a couple of options, that I know, for doing system identification with Python.
1) Using Sicpy.curve_fit() . I have a tutorial for using it with first and second order linear systems. You have to propose a system model as well as a set of initial values, a timeseries input signal (fan speed) and a timeseries output signal (temperature). Take a look of the tutorial and let me know if you have any doubt. The link.
Note. I don't recommend to use my GA-enhanced version of curve_fit() yet although you will find the link there as well.
2) Py-SINDy. This is a great package from Steve Brunton research group (University of Washington) for identification of nonlinear systems. You would have to express your linear system as a differential equation for using this tool. Here the link.
Finally, you can use this other part of my tutorial to run simulations with control systems and transfer functions. I have to admit the simulator needs to be refactored, so let me know if you cannot make it work.
Have a great day!
Saludos,
Juan.