Simscape: look ma, (almost) no equations!

I studied control engineering in school, but like many people, I found only a few occasion to directly use what I learned in school.  If I had attended this lecture live while I was at Caltech, I might have been more inspired and finished my Ph.D in control and mechanical engineering, but what is the use of thinking about "what might have been" 20 years earlier?

What little control engineering I got to do in my otherwise mostly an embedded software engineering career has been in environment control: temperature, humidity, air flow, etc.  While I still dabble in other control applications like power, motor, pneumatic, there is a universal need for this HVAC-like applications in biotech industry, and therefore pays the bills.  The control engineering method I learned in school starts with simple schematic-like diagram of the device under control, and the simplest possible relevant governing equations: for example temperature and heat flow in thermal domain, which works out to a bunch of RC-circuit like nodes.  While a single first order system is a control engineering 101 lecture material, the real world is rarely that simple: if you put a few such first order systems together, sufficiently complex behavior will result.  To even describe such system, the governing equation balloons out due to the number of states, and manually manipulating the equation before deriving a LTI (linear time invariant) system equation gets unwieldy until all the equation can be beautifully consolidated to the final state-space form: dX/dt = A X + B U; Y = C X.  But since the controller is implemented in a digital form, the 1st order differential equation is then promptly converted to the difference equation: X[k+1] = A X[k] + B U[k].  The linear algebraic formalism allowed advanced MIMO (multi-input-multi-output) controllers that can control the system from incomplete/noisy measurements.  For example, this is a block diagram of closed loop feedback regulator (tries to achieve the commanded value of the output) with an estimating observer:

A standard way to deal with noisy measurement is Kalman filter.

While great in theory, I found that complex systems (systems with say 5 states or more) are difficult to control using these LTI controllers, and made simplifications reduce the number of states before I start working with the state space representation in my ancient copy of Matlab controls toolbox, running on Windows XP.  But thanks to MSFT's finally killing the Windows XP, I just upgraded to Matlab and Simulink 2015a, and found many new things, including the Simscape--which lets you DRAW the system model (look ma, no equations!) and simulate the response.  It almost feels like cheating, but since I am not in school any more, I will explore in this blog entry how control engineering has become slightly easier since I left school.

TEC (thermoelectric cooler) Simscape model

In biotech, TEC is the temperature cooling device of choice, because the heat transfer requirement is relatively low and its non-moving part nature.  I was hoping that Simscape already includes the TEC model--if not in the fundamental models, perhaps in the SimPower or SimElectronics, but alas, it did NOT.  I ran into an article where the author built a Simscape model [Caroff, T., et al. "Transient cooling of power electronic devices using thermoelectric coolers coupled with phase change materials." 2013 19th International Workshop on Thermal Investigations of ICs and Systems (THERMINIC)] but did NOT release that Simscape model to the public domain.  So I wound up creating a custom Simscape model--not exactly a gentle introduction into Simscape.

Peltier governing equation

A TEC can be driven with either positive or negative current (with respect to the factory indicated ground terminal--usually colored black); it will move heat from one side of the TEC to the other depending on the polarity of the current.  In dealing with TEC, I found it helpful to adopt a notion that the side that gets hot when positive current is applied is the "hot side", and the heat flow direction is from positive in the cold to the hot side (yes, it's counter-intuitive, but you will see it is notation-wise simpler).

Q is a non-linear function consisting of 3 terms:

1. Peltier effect, driven by the current, which is non-symmetric across current = 0
2. Self-heating due to current flow (resistive heating), driven by the current
3. Conduction from the hot to the cold side

The governing equation for the TEC's heat moving capacity is:

Q = Sm Th Itec + 0.5 Rm Itec^2 + Km (Tc - Th) for Itec >= 0

    = Sm Tc Itec + 0.5 Rm Itec^2 + Km (Tc - Th) otherwise

Staring at the equation, I am relieved that Q is continuous at Itec = 0.  The thermal gradient is positive from Tc to Th, as I had established in the convention above.  It also makes sense that the electrical self-heating should generate heat inside the TEC.  To cool (Q < 0), Itec must be < 0, and has to overcome the resistive heating and the thermal gradient (because the Tc is usually > Th when cooling).  Clearly, the TEC efficiency will suffer drastically if the "cold side" is not adequately cooled with a heatsink.  Because the self heating grows as quadratically with the current, a TEC will eventually reach a runaway condition if it pulls a lot of current for extended period.

The Peltier effect also works the other way, and produces a "back-EMF"-like voltage across the TEC's positive and negative terminals (V+ - V-): V = Sm (Th - Tc) + Rm Itec, which can be rearranged to Itec = [V - Sm (Th - Tc)] / Rm.

A custom Simscape component for the TEC

Creating my first Simscape package

The Simscape documentation said:

Simscape files must reside in a +package directory on the MATLAB path.

After a fresh install, Matlab ~/Documents/MATLAB and put that in the path.  So I just created my own Simscape package folder there:

~/Documents/MATLAB$ mkdir +henryscape

Here, I created a new file TEC.ssc and wrote the following code:

component TEC % Thermoelectric cooler: 2.0	The first comment line after "component TEC" is the "user friendly name" to be displayed in the model view.  "2.0" after the colon
% Move heat with electrical work. % More comments not shown in the interest of space	The rest of the comment shows up as the component explanation.
nodes %conserving ports p = foundation.electrical.electrical; % +:top n = foundation.electrical.electrical; % -:bottom Th = foundation.thermal.thermal; % H:top  Tc = foundation.thermal.thermal; % C:bottom end	These ports show up in the block diagram.  The comments immediately after the domain declaration controls the displayed name and location of the ports in the block diagram
parameters S = {0.08, 'V/K'};  % Siebeck constant K = {2, 'J/(K*s)'}; % thermal conductivity; should use [W]? R = {2.2, 'Ohm'};   % electrical resistance end	Parametrizing the component's property enhances reuse.  In my TEC component's case, the 3 parameters were captured in the governing equation given above.  The default values are for the TEC I happen to be using, but easily changeable when I double-click on the component in the model editor.  I hesitated a bit when picking the unit for the thermal conductivity: Watts per Kelvin [W/K] would have been more concise, but I stuck with [J/(s*K)] because I saw the DC motor component example use [J/s] when it could have used [W].
variables(Access=private) i = { 0, 'A' };  % Current v = { 0, 'V' };  % Voltage qK = { 0, 'J/s' };% conductive heat flux  qR = { 0, 'J/s' };% resistive heat generation  qS = { 0, 'J/s' };% Siebeck heat flux end	I actually don't know yet the reason for choosing different access level; just copied from a Simscape component example.  Could I have used 'W' instead of 'J/s'?  Yes.
function setup if S <= 0 pm_error('simscape:GreaterThanZero', 'S') end if K <= 0 pm_error('simscape:GreaterThanZero', 'K') end if R <= 0 pm_error('simscape:GreaterThanZero', 'R') end end	Parameter validation.
branches i : p.i -> n.i; qK : Th.Q -> Tc.Q;  qR : * -> Th.Q;  qR : *-> Tc.Q;  qS : Th.Q -> Tc.Q; end	qK (conduction always goes from hot to cold).  Keep the sign convention this way and also keep your sanity. Resistive heat is generated from inside the TEC and flows outside. Peltier heat flows from the hot side to the cold. branches specifies the relationship between the Through variables
equations v == p.v - n.v; v == S * (Th.T - Tc.T) + R * i; qK == K * (Th.T - Tc.T);  qR == 0.5 * R * i * i; if ge(Th.T, Tc.T) qS == S * Th.T * i; else qS == S * Tc.T * i; end end % equations	Equation specifies the relationship between the Across variables. qR is always NEGATIVE so that heat always flow from the TEC to the outside.  But this is OPPOSITE of the above definition. qS is based on the HOTTER side temp.  It is NON-linear. "let" keyword would allow me to describe TEC heat flow with intermediate terms for comprehension.  end keyword is required at the end of a let-in-end statement.
end %component	I miss the curly braces...

Compiling my custom package

I also wanted a TEC icon for the block, so I found a TEC image (below) on the web and put it right next to TEC.ssc.

I like this picture because it shows the heat "going up" as I've written the equation.  To actually generate the package, I ran this command in the Matlab shell:

>> cd ~/Documents/MATLAB
>> ssc_build henryscape
Generating Simulink library 'henryscape_lib' in the current directory '/home/henry/Documents/MATLAB' ...

After this, I find a file henryscape_lib.slx in ~/Documents/MATLAB folder, and when I double-clicked on the lib file in Matlab's folder, I saw the component appear in the Simulink libarry viewer:

Note that the + and the H nodes are on top, and the - and C nodes are on the bottom, as intended in the source code.

Using the TEC component in a Simscape model

I create a new Simscape model with this Matlab command:

>> ssc_new

This populates a new model with components that every Simscape model needs, as shown below:

Note the model is using the recommended solver for Simscape: ode23t.

I then dragged the TEC component from the library viewer shown in the previous section into this model, and then added an ideal current source controlled by a square wave, and measured the temperature on the hot side of the TEC, using this model:

The resulting plot of the Th - Tc (the temperature difference between the TEC hot and the cold side) is shown here, overlaid against the +/- 1A current being driven from the current source:

Even with an ideal current source, the Peltier effect prevents the current from changing discontinuously.

Next step

I will connect this TEC component to the part that I am trying to heat/cool, and also to the heatsink that sits under the TEC.  I saw a Simscape tutorial that fits parametrized model to openloop data.  In that scenario, some of the parameters are left in symbolic form that can be changed from Matlab script.  For example, a thermal capacitance (mass) of the part I am trying to cool is called wallMass in the thermal mass block property, as shown below:

Ultimately, I will throw in a real-time controller that acts on the measured temperature of the part that I am trying to heat/cool.