My only outdoor activity, and the majority of my daily exercise is biking to and from work. I used to wear an ANT+ Garmin chest strap paired with the ANT+ iPhone 4 adapter to measure my heart rate and estimate my caloric expenditure. The readings I got in the Digifit report was accurate when it worked, but had a few problems:
Initial sync with the my iPhone took up to a minute.
Unless I wore my iPhone on my arm to keep it close to the sensor, it was prone to losing connection. The monitoring application I used--digifit--was not smart enough to recover from intermittent connection drop, so I would sometimes lose more than an hour of riding record.
I got fed up, and stopped using it.
Recently, I started working at Jawbone--specifically on the wearable band firmware. One of the perks is a free UP3 band, which monitors sleep, estimates steps, calory, and resting heart rate.
Because it is always running, I don't have the hassle of maintaining connection with my cell phone (except for a daily sync with the cell phone). I think it is a wonderful hardware, but it is not sensitive enough to accurately monitor my biking activities, because my hands are moving very little when I am on my bike.
One workaround suggested to me is to put it in my pocket. Since my thigh is moving up and down when I pedal, it can count steps while in my pocket. Another workaround was to get an UP MOVE, and clip it to my bike shoes.
But in truth, I don't like wearing anything on my wrist, so I thought: why not just wear the band on my ankle? Then it can detect my heart rate just like when it is on my wrist! The problem is that the band is too short to wrap around my ankle, which is 8 3/4 inch, as measured with a string around my ankle. If my ankle circumference had been half inch shorter, a fully extended UP3 (see the picture below) would probably wrap comfortably around my ankle. I tried it, and found it too tight. When the clasp is at the longest position, there is actually more than a half inch left from the end of the rubber band--I just could not get there because of the bio-impedance sensors in the way, as shown below:
WARNING: I do NOT speak for Jawbone in describing my hacking below. Doing what I did will probably void your warranty, so please don't blame me or Jawbone if you break your UP3 following my example. That said, UP2 doesn't have a bio-impedance circuit in the rubbery band, so you probably will not damage it.
So I took off the adjustment clasp with a pair of pliers, and flipped it around 180 degrees to have it FACING the other end of the clasp pair, as shown below:
Further bending of the metal made it feel comfortable around my ankle, as you can see below:
I am happy with the result: the UP3 band is working as well as it did when it was on my wrist, and of course now accurately counts my pedaling. I was lucky: if my ankle circumference had been even half inch longer, the band would have felt too tight.
3 month update
The clasp is looser than before, and regularly comes off during sleep. I don't recommend this hack for others. I took a plier and tightened up the clasp a little, and now it's very secure.
UNFINISHED WORK published only for collaboration with others
I was recently introduced to compressed sensing, so I am going to teach myself how to use it to my problems. Instead of struggling with CS math treaties, I will go through as many programming examples as possible, like this one aired on NPR, which makes a point that piano music can be sparsely coded and decoded with a basis consisting of piano notes, but not (very well) human speech. I liked this tutorial because it taught me about taking advantage of duality when designing the sampling strategy: if the signal is sparse in time, sample randomly in frequency; conversely, if the signal is sparse in frequency space, sample randomly in time.
Introduction
Imaging in next gen sequencing
The next-gen DNA sequencing derives its high thoughput by observing as many randomly fragmented DNA strands in a given reaction area as possible. But since DNA strand is a molecule (polymer, more precisely), it is damned near impossible to see (which is where the 3rd gen sequencing comes in), a DNA strand is amplified to a cluster big enough to be seen with high resolution collection optics, as shown in the cartoon schematic of a Solexa/Illumina cluster:
If all goes well, a high resolution camera can observe multiple clusters fluorescing with a given dye. By aligning the images for different dyes, a base calling algorithm can classify a given cluster to be fluorescing with one of the possible dyes at a given cycle, as shown in this example, where we can manually call the bases for 2 clusters (circled):
Of course, the reality is never this clean. To observe as many clusters in 1 frame as possible, we want to size the cluster to span as few pixels as we can get away with, as in this example of a portion of "T dye" image of an early Solexa run rendered in ImageJ:
An average cluster's PSF in these runs is on the order of a micron, so a pixel is on the order of a 200 nm.
Problem: high resolution imaging is expensive and difficult, but low resolution imaging cannot do the job
Focusing and aligning optics system down to an order of 100 nm requires many things to work exactly as they should. Suppose that I have an imaging system that can resolve 1 um (which is not easy either!). Even if I could focus the clusters to this poorer resolution camera (I can't), the clusters would be so smeared (binned) that I cannot pick up a putative cluster from the resulting image like this:
I want to see if the original image is sparse enough to benefit from compressed sensing.
Simulating random imaging in Matlab
Single pixel camera experiment showed that a sparse (in space) image can be reconstructed from randomly oversampled single-photo-detector data. The random sampling was done with a DLP, which had sufficient spatial resolution. In Matlab, taking ONE random sample means to convolve the signals with a randomly chosen basis vector--and DFT (discrete Fourier transform) is as good a basis as any.
1D DFT
To make sure I understand the operation, let's fall back to the definition of a DFT given a vector x[k], k = 0 through N-1:
X[n] = SUM_k { exp(2i n k/N) x[k] }, k and n = 0~N-1
= [exp(2i n 0/N) ... exp(2i n (N-1)/N)] * [x[0] ... x[N-1]]'
So X can be computed with a matrix and vector multiplication: X = A x
[ X[0] ] [... ... exp(2i 0 (N-1)/N)] [ x[0] ]
[ X[n] ] = [... exp(2i n k/N) exp(2i n (N-1)/N)] * [ ... ]
The square matrix A matrix above is the DFT matrix, obtained in Matlab with the function dftmtx(N). To invert the transform, complex conjugate of the matrix is required, because the definition of the inverse transform is identical to the forward transform, except for the negative sign in the exp argument.
x = conj(dftmtx(N)) * X .* (1/N)
An image is commonly represented as a 2D matrix, but is in fact commonly stored as 1D vector in memory. So 1D can be used without loss of generality. BUT, it is difficult for people to "see" an image unrolled to 1D, so a 2D treatment is useful.
Random sampling in 1D frequency space
In the Matlab CS tutorial mentioned above, random sampling in discrete frequency space is accomplished by randomly selecting K rows of the DFT matrix (which is square), which is extremely easy to do in Matlab:
qN = randperm(N); randN = dftN(qN(1:K), :);
y = randN * truth + v;
where "truth" is the unobservable true state that we are trying to estimate, and v is an additive noise term that is bound to be present all real world measurements--usually set to 0 at first pass through the theory. Actually implementing such sensing hardware is extremely difficult: not the noise injection part--which the nature does an excellent job of--but the construction of a hardware that will pseudo-randomly mix the spectral components of the true states into under-sampled measurements.
The result of the round trip is the maximum value of the difference between the original image and the round trip image, which turns out to be 0.
Flattening an image to 1D vector
For compressive sensing, 2D data is unwieldy, because the literature and examples all use 1D data. A silver lining is that in computer memory, multi-dimensional memory is already flattened out to a vector, and flattening a matrix to a vector in Matlab is trivial.
img_flat = img_orig(:); % Flatten an image into a vector
BUT flattening a 2D DFT to a 1D will NOT yield the same spectra as the 1D transform of the flattened image: fft2D = fft2(img_orig); fft2D_flat = fft2D(:); fft1 = fft(img_flat); >> max(fft2D_flat - fft1) ans = -8.7022e+07 + 2.7322e+07i
Sparse random sampling of image
The example image I am using is 1024x1024 pixels, with each pixel representing roughly 200 nm squared surface area. Suppose I get a sensor with 1.6 um resolution (8 times worse than the original).
I've used the full-state Luenberger observer successfully in my control problems, using Digital Control of Dynamic Systems, 2nd edition, Franklin et. al., chapter 6 as my reference. The key block diagram that allows for control commands and current estimation is this:
Although this structure does NOT have error integral portion and may therefore have a non-zero steady-state error, the later part of the above mentioned chapter does give the modified solution, and besides in my experience, integral control ALWAYS requires anti-windup correction, so the simpler feedback gain formulation given above is a good first attempt at a non-full-state-feedback MIMO (multi-input, multi-output) control problem.
I in fact learned about the observers from Control System Design Guide, George Ellis, who to his credit, tried to make a dent against the perversely pedantic approach to control engineering, by creating an excellent control system simulation software (for nearly free!). In my opinion, for SISO (single-input-single-output) design, the Visual Model Q is superior to Matlab/Simulink because it helps (forces?) me to understand what is going on under the hood.
Unfortunately, I have to solve MIMO control, estimation, and optimization problems, and Matlab toolboxes are a faster way to tackle those problems. And after avoiding Simulink all these years, I finally bought into it recently. Part of the entry fee into the Simulink was implementing the above block diagram in Simulink. Here is my first attempt (where problem specific--proprietary--labels are erased):
The predictor that is driven by U (control input into the openloop model) and Y (measurement, or observation) is fairly standard, with the predictor's system matrix A modified to (A - Lp C), and the B matrix B modified to [B Lp], where Lp is the predictor feedback gain. The 2 inputs on the left is the R vector to follow, the 2 independent commands in this problem.
But there are some important differences:
I did not put in Nx, because in this problem, Nx = I. But adding a matrix gain between R and the difference junction (circle above) is easy; its gain parameter will just be the Nx matrix, available in Matlab workspace.
Feedback is calculated only on the subset of the full Xest, so I use another gain (a 2x5 matrix consisting of only 2 rows from a 5x5 identity matrix).
Because controller saturation is so important, I explicitly modeled it. In fact, when the saturation limit varies as a function of the states, I modeled that relationship as well, and feed Simulink's dynamic saturation block.
Future work
Note that I don't deal with the system or the measurement noise, by using the Kalman filter, for example. The full-state predictor/estimator in fact averages out noise nicely, but the Kalman filter will dynamically adjust the optimal weight for the current vs. past estimates.
Whenever I come across an Xilinx FPGA tutorial that says: "make sure your design is free of validation errors/warnings", I get skeptical about the writer's experience level, because I've NEVER come across a non-trivial Xilinx design without hundreds of warnings, and a dozen or more critical warnings. The best I can hope for is to be free of outright errors. Consider the following critical error in my Vivado design:
[BD 41-967] AXI interface pin /axi_i2s_adi/DMA_ACK_RX is not associated to any clock pin. It may not work correctly.
...
I've been living with this warning ever since I got a Zedboard. Recently, I finally looked it up on the Xilinx forum. These are false alarms, according to a Xilinx forum discussion.
In previous blog entries, I created custom Simscape components for TEC and forced air heatsink for a thermal control system. I built an open-loop system model with a current source, and now I want to find reasonable values for the model's parameters using openloop data. While I've done this exercise many times in the past, I always did it within Matlab itself rather than in Simulink. In an effort to catch up to the modern Matlab workflow, I will document my new work procedure.
Create input/output ports for the model
I already have my model (creatively called "take1") open, but it is currently missing output signal. I dragged an out1 port from the Simulink Library Browser --> Simulink --> Sinks, and changed the port property to show port name and number. Similarly, I dragged an in1 port from the Library Browser and connected to the ideal current source in my model.
Bring the openloop data into Matlab workspace
My data is already in a CSV file in my Matlab folder (~/Documents/MATLAB), so I can easily bring it into the workspace by right clicking on the file in the Current Folder view --> Import Data... --> Import. Now my Matlab Workspace is choke-full of column vectors.
I believe the csvread() command would work just as well. In there, I identify the vectors I want to use for this experiment:
input1: TECValue, commanded value to the current generating FW
output1: CenterTemp
ouput2: I have several thermistor readings, which I should average, like this:
>> center = mean([Sensor5 Sensor6 Sensor9 Sensor10], 2); >> center = center(4215:6000); % Just deal with relevant data
output3: HeatsinkTemp
The data was captured at 2 Hz, but the time was not recorded. So I just make it up with this Matlab command:
>> time = (0 : 0.5 : 0.5*(size(center,1)-1))';
Note the transpose (') operator to keep the time as a column vector, just like the other openloop data columns. I suppose I could have defined a sampling frequency as a constant, but let's just keep moving for now.
Create an experiment
In the Simulink model window, select Analysis > Parameter Estimation. This action opens a new session, Parameter Estimation - take1, in the Parameter Estimation tool. Click "New Experiment" in the toolbar to bring up "Edit Experiment: Exp" wizard--"Exp" is the default name of the experiment.
Select all output ports I created in the model, and supply the openloop data
output 1: [time, CenterTemp(4215:6000)]
output 2: [time, center]
output 3: [time, HeatsinkTemp] I ignored the heatsink temperature because it is actively controlled.
In the "Inputs" category, select only some of the rows: [time, TECValue(4215:6000)]
Select Initial States. Optional?
Select parameters: Simulink found all parameters in my model: those that were left in symbolic form. I checked all of them. If necessary (usually a good idea), the parameters should be constrained to sensible bounds. when expanded, the parameter tab shows the minimum and maximum bounds. For almost all physical parameters, the minimum should be some small value.
Clean up data
To view the imported data, click Add Plot on the Parameter Estimation tab and select the experiment name--"Exp" here.
Since the open loop data covers many different cases, here is a way to restrict the data range: in the "EXPERIMENTAL PLOT" tab, click "Extract Data" button, and specify the "Start/End Time", as shown below:
Then I clicked "Save As", which created another experiment "Exp1". I renamed Exp1 to "Exp TEC". In the EXPERIMENT PLOT tab above, there are other ways to clean up the experiment data:
Remove outlier, offset
Scale, filter data
Estimate parameters
When I click on the "Estimate" button shown below, the GUI runs optimization loop, and adjusts the parameters until the score (sum of the squared error) stops converging.
The EsitmatedParams "Scales Value" shows that my initial guess for one of the parameters was particularly poor. Going back to the Experiments browser --> Edit --> Parameters, the parameter values at the end of the estimation iteration are shown, and I can just read them off.
Clicking on "Plot Model Response" button above overlays the simulation result against the experimental data.
The fit is not all that great, suggesting that my simple model is missing at least 1 first order dynamics between the TEC and the out1, and maybe the TEC data I have is somewhat off. Rather than complicate the model further, I let the parameter estimator vary the TEC parameters to obtain a better fit, as shown below:
If I had more openloop data, I could validate the model + parameter against it.
Next step
I have to design a feedback and feed-forward controller for this system. The feedback controller will be relatively easy, but the feed-forward tends to have lots of heuristics and tuning knobs.
In a previous blog post, I simulated a temperature controlled device that uses TEC to move heat in either direction. Another common practice in thermal control is to thermally glue a heatsink on the "cold side" of the TEC. When cooling an area thermally glued to the "hot side" of the TEC, the "cold side" gets hot, and reduces the cooling efficiency; if the "cold side" continues to get hotter, the TEC will go into a runaway condition. The solution is to increase the effective surface area of the TEC cold side with a heatsink. But as you see in your own laptop, natural convection occurring on the larger surface area alone cannot move the many Watts the modern TEC can move and need the help of forced convection cooling; air cooling is cheaper than liquid cooling.
Forced convection constant is a non-linear function of the fan speed
Regardless of whether forced or unforced, all heatsinks have thermal mass (capacitance). So an easy way to think about the problem at a high level is to consider the change in the thermal resistance from the heatsink fins to the ambient (air if air cooled, or liquid if liquid cooled). If I plot Tsink - Tambient (the temperature of the heatsink relative to the ambient) with DC fan at max speed, I may see a classical RC circuit cooling behavior, like this:
The expotential decay portion is the exp(-t/RC) curve, so I can plot ln(Tsink - Tambient) and read off the negative slope to estimate the 1/RC. Intuitively, higher DC fan speed yields higher 1/RC--faster cooling, BUT 1/RC does NOT scale linearly with the fan speed.
Heatsink convection vs. air speed
A simplified model of convection is Q = h (Tsink - Tambient), and the problem boils down to estimating the convection constant for various conditions. I remember the scary integrals in the heat transfer class, which are often done in FEA these days. I didn't do well in fluid mechanics or heat transfer, but was able to follow along a heat transfer textbook: Engineering Heat Transfer 2nd Edition, by Rathore and Kapuno. On page 563, there is a heatsink fan configuration that looks similar to mine (array of cylinderical rods facing the air), so I started there. For forced air convection, the h term is Nu ∙ kf / D, where
kf is the thermal conduction coefficient of air film around heat fin. Typically, we would just look up the value at an average temperature of (Tsink + Tambient)/2. Between 26.3 ~ 30E-3 W/m∙K for 300~350 K.
D is the rod diameter, say 2 mm.
Nu is the Nusselt number that depends on Reynolds number (which in turn depends on the forced air speed). For my configuration, the textbook says it is 0.27 ∙ (u_air D / ν)^0.63 ∙ Pr^0.36 ∙ (Pr∞ / Prsink)^1/4, where
u_air is the bluk air speed through the heatsink,
ν is the air viscosity, ranging 15.89 to 20.92 m2/s for air between 300 to 350 K,
Pr is the Prandtl number: 0.707 at 300 K and 0.700 at 350 K.
The Nu therefore works out to something like 0.7E-3 u_air^0.63
Note that Nu depends the most on the air speed uair, namely, non-linearly as a power of 0.63. With this Nu, you are supposed to solve for the temperature of the air leaving the heatsink fins: Texit = Tsink – (Tsink – Tambient) exp[-π D N h / (ρ u_air N_T p_T C_p)], where:
N is the total number of fins
h was calculated above,
N_T is the number of rods in the transverse plane
ρ is air density, about 1 Kg/m3 around 300~350 K
p_T is the fin spacing (pitch)
C_p is the air's heat capacity (about 1 J/g∙K)
After substitution, you arrive at Texit = Tsink – (Tsink – Tambient) exp[Ca / u_air0.37], where C_a is depends on your situation. What this equation says is that the exit temperature at the heatsink is almost the Tambient at all air speeds except the really slow case, in which case the exit temperature is almost the Tsink, driven by this example of the exponential term:
Intuitively, convection from a heatsink is NOT from Tsink, but rather a temperature that is closer to Tambient. How close depends on the heatsink shape and air speed.
Finally, what you really want to know is the aggregate heat flow from the heatsink to the ambient air. For the N number of rod case, Qconvect = N (h π D L) ∆T = π N ∙ Nu ∙ kf ∙ L ∙ ∆T, where:
∆T is calculated above--practically speaking something like 90% of (Tsink - Tambient),
L is the rod height
Since Nu is proportional to (u_air D / ν)^0.63, Qconvect ~ some_number ∙ u_air^0.63 ∙ (Tsink – Tambient) for all but the very low air speeds.
Now that you understand, keep the simplest 1 equation
If the textbook was completely applicable for my case, I should be able to predict the 1/RC slope for all forced air cooling based on the airspeed--which according to the DC fan manufacturer scales linearly with the fan speed. But instead of u_air^0.63, I found a factor of u_air^0.8 to be a better fit. So instead of losing sleep over the discrepancy, I will just use the power of 1/0.8 = 1.25 when calculating the desired fan speed. That is, my forced air heatsink is simply characterized by
where 0 <= duty <= 1 ; that is, consisting of the forced and natural convection terms.
Custom Simscape component for the forced air heatsink
I explained the steps to creating a custom Simscape component in a previous blog. Here is the Simscape code for the above heat flux equation:
component Heatsink
% Heatsink: 2
% If the duty cycle is s, where 0 <= s <= 1,
% q = Csink (h100 * s^0.8 + h0) (H - C)
%
% where the following parameters are treated as constants for this simple
% model:
%
% * h100: convection constant at max fan duty cycle [W/K]
%
% * h0: natural convection constant [W/K]
nodes %conserving ports
H = foundation.thermal.thermal; % H:top
C = foundation.thermal.thermal; % C:bottom
end
parameters
h100 = {0.032, 'W/K'};
h0 = {0.0019, 'W/K'};
end
inputs
s = {0, '1'};% s:bottom
end
variables(Access=private)
q = { 0, 'W' }; % heat flux
end
function setup
if h0 <= 0
pm_error('simscape:GreaterThanZero', 'h0')
end
if h100 < 0
pm_error('simscape:GreaterThanZero', 'h100')
end
end
branches % Relationship between the Through variables
q : H.Q -> C.Q;%In heatsink, heat flows hot to cold
end
equations % Relationship between the Across variables
q == (h100 * s^0.8 + h0) * (H.T - C.T);
end % equation
end %component
And while I may yet get tired of such frivolities, I got a picture of a heatsink from Bing Images and put it right next to the source code during the compilation of my simscape library, so that this custom component may show up as a picture.
Simulate cooling a hot block with the heatsink
As a sanity test, set the initial temperature of the heatsink at 100 C, when turn the fan on at 50% duty cycle. Here is the Simscape block diagram for it:
And here is the 10000 second simulation result:
As you can see, this heatsink is not all that powerful.
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:
Peltier effect, driven by the current, which is non-symmetric across current = 0
Self-heating due to current flow (resistive heating), driven by the current
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.
i n k/N) x[k] }, k and n = 0~N-1
