Thermoelectric Modules in Mechatronic Systems: Temperature-dependent Modeling and Control

Active thermal control is crucial in achieving the required accuracy and throughput in many industrial applications, e.g., in the medical, high-power lighting, and semiconductor industry. Thermoelectric Modules (TEMs) can be used to both heat and cool, alleviating some of the challenges associated with traditional electric heater based control. However, the dynamic behavior of these modules is non-a ne in their inputs and state, complicating their implementation. To facilitate advanced control approaches a high fidelity model is required. In this paper an approach is presented that increases the modeling accuracy by incorporating temperature-dependent parameters. Using an experimental identification procedure, the parameters are estimated under di↵erent operating conditions. The resulting model is used in a feedback linearization approach to linearize the system, facilitating the use of advanced linear control techniques. Moreover, it presents an observer based approach to reconstruct state information in cases where sensor placement is limited. The resulting framework forms a complete approach to temperature-dependent modeling and control of thermoelectric elements.


Introduction
Advanced thermal control is crucial to achieve the required specifications on accuracy and throughput in many industries, especially in the medical, high-power lighting, and semiconductor industry. For example, in the medical field, diagnostic platforms are used to process extremely small fluid volumes, e.g., blood or saliva samples [1]. Hand-held devices are designed to significantly reduce analysis time and reagent costs. The temperature of the fluid volumes needs to be accurately controlled during a wide range temperature cycle. In particular, in PCR amplification, which is presently used to test for active infections in patients during the COVID-19 pandemic [2], it is crit-ical to quickly cycle the sample over a wide temperature range. Another example is in high-power LED lighting, where LEDs o↵er superior e ciency compared to traditional lighting. However, due to their limited volume the heat that is produced is of significant intensity and must be actively controlled to achieve su cient light quality and an increased lifespan [3,4]. Finally, in the semiconductor industry, wafer scanners are used to produce integrated circuits that need to achieve a positioning accuracy of nanometers, where the current performance of these high precision systems is often limited by thermally induced deformations. Therefore, thermal control is an important aspect in their mechatronic design [5]. In [6,7], the benefits of active thermal control are already confirmed for this class of systems.
TEMs are not limited to heating, o↵ering major advantages over resistive heaters that are often used in thermal control. TEMs are compact, have no moving parts and facilitate both heating and cooling. This lends them particularly well to appli-cation in compact systems, e.g., small scale PCR amplification devices [8]. In mechatronic systems they can be employed to locally provide both heating and cooling.
The thermodynamics of a TEM are non-a ne, containing both state dependencies and input non-linearities, complicating the controller design. Standard linear control methods (e.g. PID control) could be unreliable because stability, robustness, and performance cannot be guaranteed. Therefore, nonlinear control methods can be considered. In recent literature several methods have been developed. In [9], a linear-parameter-varying approach is used to control the nonlinear system, which linearizes the nonlinear system at di↵erent operating points. For each operating point a di↵erent controller is synthesized. In [10], a sliding-mode controller is used, which applies a state feedback. The state feedback ensures that all trajectories move towards a stable sliding manifold. In contrast, in [5], the nonlinear system is partially linearized using feedback linearization by creating a new input that has linear input-to-output (IO) dynamics. This facilitates the use of conventional linear control approaches. However, this approach does not include temperature dependency in the model and linearization controller. This constrains the applicability to situations were the temperature remains fairly constant, limiting the usefulness for the industrial applications considered in this paper where a large temperature range is desired.
In related work [11], a linearization approach is used to control the cold side of a TEM for a large temperature range from 5 to 80 C. However, it is shown that the feedback linearization yields some residual nonlinear dynamics, that might cause instabilities. Moreover, in both [5] and [11] an observer design is recommended since it is often impractical or infeasible to install temperature sensors around the point of interest (POI). For example, in a diagnostic platform the fluid temperature must be accurately controlled but placing a temperature sensor in the fluid is undesired due to hygiene constraints. In [12,13], results using an Extended Kalman Filter (EKF) are presented and the stability proof for this class of systems is involved [14].
In [15,16] it is shown that the physical parameters of Thermoelectric modules vary with temperature. However, in previous work [12,5,11], the model parameters are typically considered as constant. In sharp contrast, the main idea pursued in this paper is that of temperature-dependent mod-eling and control. This is achieved by including temperature-dependent parameters in the model, feedback linearization, and observer design.
This paper provides a complete framework for temperature-dependent modeling and control of thermoelectric elements. A dedicated experimental identification setup is used to obtain a highfidelity temperature-dependent model, presented in earlier work in [17]. By leveraging the temperaturedependent parameters in the model, the system is linearized using feedback linearization that is valid over a wide temperature range. And to facilitate situations where sensor placement is limited, an observer based approach is presented that is capable of providing an accurate state reconstruction using a limited number of temperature sensors. The main result of this paper is a unified temperature-dependent modeling and control approach for thermoelectric elements, with the following detailed contributions. C1 A modeling approach to obtain an accurate temperature-dependent model for a TEM.
C2 A systematic estimation procedure for temperature-dependent parameters over a wide temperature range.
C3 Incorporating temperature-dependent parameters into a feedback linearization approach, with a stability guarantee.
C4 Incorporating temperature-dependent parameters into an observer design. Accompanied by a Lyapunov-based stability proof.
C5 Experimental validation of the presented approach, from modeling to control.

Problem formulation
In this section, a motivation and problem formulation for modeling and control of thermoelectric elements is presented. Subsequently, the experimental setup that will serve as an illustrative benchmark setup throughout this paper is presented.
In view of control and to facilitate the implementation of both accurate linearization methods and observer design, a high-fidelity model of the TEM is required. In literature, often a limited operating temperature for the TEM is considered, allowing the model to be simplified by using temperature independent parameters. In this paper, in particular to facilitate the application of peltiers in the medical industry [8,18], a significantly larger operating temperature range is considered, e.g., from 5 to 80 C, necessitating the inclusion of temperature dependency [19,21] in the simulation model.
Moreover, it is important to take the increased temperature range into account during the control design. Previous results have shown that a feedback linearization approach is promising but can lead to inaccurate results when neglecting temperature dependencies. Concurrently, it is important to provide a suitable stability proof for the entire operating range. Finally, to facilitate situations were sensor placement is limited, a suitable observer based approach must be constructed.
Summarizing these aspects leads to the following problem formulation. Provide a framework for temperature-dependent modeling and control of thermoelectric elements that facilitates their application in a wide range of operating conditions.

Experimental setup
In this section, the experimental setup that is used throughout this paper is presented. The methods and techniques presented in each section of this paper are validated using this experimental setup.
The setup, shown in Figure 1, consists of two thermoelectric elements clamped between a steel bottom plate and an aluminum top plate. Both elements are cooled by convective cooling using an aluminum heat sink attached to the top plate. The experimental setup incorporates 2 TEMs to provide additional challenges by introducing interaction. More precisely it uses TEM 2 en 3 from Section 4 since they are the most similar. The setup is constructed such that it presents a su ciently challenging benchmark for modeling and control and will be referred to as the application setup throughout.
To measure the temperature, voltages and current in the TEM setup a CompactDAQ by National Instruments is used. To facilitate temperature measures using the thermistors with Negative Temperature Coe cient (NTC), a Wheatstone bridge is used that converts the resistance measurements, and thereby the temperature, to an electrical potential. Moreover, since the identification method proposed in Section 4.1 relies heavily on accurate knowledge of the input current, a precision power resistor is placed in series with the TEM. The resistor is selected such that its resistance R remains constant for the operating currents. By measuring

Modeling
In this section the first principle relations describing the major thermoelectric e↵ects within a thermoelectric module are presented. Following this, a full lumped-mass model of the experimental setup presented in Subsection 2.1 is provided that is used throughout as an illustrative example.

First principles
The thermoelectric module in this paper is modeled through lumped-capacitance discretization, i.e., the module is subdivided into lumps of uniform temperature. This is illustrated in Figure 2, where the module is divided into a hot and cold side where each ceramic plate is a single lump. Conversely, some literature models a TEM on the pellet level [22], which would be overly complicated for the application considered in this paper. Moreover, it is shown that a lumped-mass discretization can yield su cient model granularity.   The thermoelectric dynamics of a TEM are described by including 3 phenomena: 1) the Fourier e↵ect Q f , 2) Joule heating Q j and 3) the Peltier e↵ect Q p . The related Thomson e↵ect [20] is considered outside the scope of the current paper.
Fourier e↵ect. The Fourier e↵ect Q f describes the energy transfer through conduction between the 2 sides of the TEM and it is given by for conduction from temperature T 1 to T 2 with K m the conductivity of TEM in W/m · K, A the area in m 2 perpendicular to the heat flow and d in m the distance of the heat flow path.
Joule heating. Joule heating Q j occurs when an electrical current flows through a resistive element, in this case the TEM, and is given by where R m is the electrical resistance in ⌦ of a single TEM and I is the current in A. It describes the energy loss of electrons moving through an electrical resistance, and it is typically unavoidable without employing superconductivity at cryogenic temperatures.
Peltier e↵ect. The Peltier e↵ect describes the occurrence of a heat flow over a semi-conductor in the presence of an electrical potential di↵erence and resulting current. Conversely, the Seebeck e↵ect describes the occurrence of an electrical potential over a semi-conductor in the presence of a temperature gradient. While they are manifestations of the same physical phenomenon, the former is described as where S m is the Seebeck coe cient of the TEM and T is the temperature at the cold/hot side. It is assumed that the Joule heating Q j , that is generated in the semi-conductors, see Figure 2b, is divided equally over the hot and cold side of the TEM since the module is symmetric. The heat balance equations for the hot and cold side are then given by where Q env!c,h accounts for any thermal interaction with the environment, indicated by the superscript env , i.e., the ambient air or neighboring lumps.

Modeling the experimental setup
To construct a finite order model that is tractable for parameter identification and control, the setup, as shown in Figure 1, is discretized into 18 lumps. Lumped mass modeling assumes the lump temperature to be uniform in each discrete lump, and a measure for the accuracy of this assumption is the Biot number that is a ratio of the internal and external resistance which should be smaller than Bi < 0.1 to achieve uniformity within the lump. Consequently, a discretization as shown in Figure 3 is obtained.
To obtain lump temperature measurements, both thermocouples (TC) and negative temperature coe cient thermistors (NTC) temperature sensors are attached to the setup. Air temperature x 9 x 11 x 10 x 12 x 16 x 17 x 15 Qc(x10, x12, u2) Figure 3: Illustration of the lumped-mass discretization of the experimental setup resulting in a network of connected nodes. The illustration also indicates the sensor placement, which is assumed to measure the true lump temperature.
By constructing an energy balance equation for each lump, an accurate model is constructed including the TEMs and any connecting elements. This results in a state-space model, where the states x = T (1,...,nx) , with n x the number of states, represent the temperature of the lumps with corresponding state equations where E n = m n c n is the thermal capacitance of the lump n with m n the mass in kg and c n the specific heat capacity in J/kgK and P Q n the heat balance equation, akin to (4), for a specific lump n. By collecting these di↵erential equations the statespace model of a system is given by is a nonlinear function depending on x and the input current u = I, C 2 R ns⇥nx is the output matrix with n s the number of temperature sensors, x(0) n 2 R nx⇥1 is the initial condition of the states. By incorporating the nonlinear input, e.g., the joule heating Q j and state-dependent dynamics in F NL (x, u) a full system model is constructed. This nonlinear input function F NL (x, u) is given, for the illustrative setup, by , with n M the number of TEMs in the setup, are respectively the thermal conductivity, electrical resistance and Seebeck coe cient of the corresponding TEM. Observe that the first term in (9) is actually a linear contribution, assuming that K m is not state-dependent. Including this linear term in (9) allows the structuring of the linear dynamics matrix A 2 R nx⇥nx as where A H 2 R n H ⇥n H and A C 2 R n C ⇥n C are considered the hot and cold sides dynamics respectively with n H ,n C the number of lumps in their respective parts of the model, i.e., n H = 10 and n C = 8 for the illustrative setup. Consequently, the states are structured as are the states associated with the thermal node network on the hot side and cold side respectively, i.e., x 1,...,8 for the hot side and x 13,...,18 for the cold side. The disturbance input matrix B 2 R nx⇥1 is also structured as where B H 2 R n H n M ⇥1 and B c 2 R n C n M ⇥1 since it is assumed that the TEMs are not coupled directly to the ambient air w, see Figure 3. The output matrix C 2 R ns⇥nx is structured as where such that the appropriate lumps are measured according to the sensor placement shown in Figure 3. Moreover, the lumps x 9,...,12 consist of the TEM ceramic plates and their state, i.e., temperature, is not directly measurable.
Remark 1. Note that the derivation of (8) is done in continuous time, a discrete time model is obtained by applying an appropriate discretization method, e.g., Euler discretization. For completeness and to allow reproducibility of the results, a full model of (8) is provided in Appendix A.

Temperature-dependent modeling
Employing constant parameters in the model (8) can yield su ciently accurate results for systems that operate in a limited temperature range, as demonstrated in [5]. For the systems considered in this paper, e.g., a blood diagnostic device that cycles over a large temperature range, this is often not su cient and temperature dependencies must be taken into account. In this section, emphasis is placed on temperature-dependent modeling, and it is shown that including this dependency can increase modeling accuracy for a wide temperature range, providing Contribution C1. Moreover, it is shown that there exists a systematic approach to obtain these temperature-dependent parameters through dedicated experiments, providing Contribution C2.
Including the temperature dependency in (8) is done by modeling the parameters S m and R m as a function of the average temperature T avg of the TEM, i.e., S m (T avg ), R m (T avg ), where Remark 2. While the conductivity of the TEM is also considered temperature-dependent in some literature [23] in this paper this was not deemed beneficial within the current temperature range and it is taken as a constant value.

Identifying temperature-dependent parameters
To include the temperature dependency of R m and S m requires an identification procedure that can determine the value of these parameters over the required temperature range. This is done by measuring the electrical potential V M required to induce a fixed current I M in a single TEM. The total voltage is given by where I M is the current output in ampere of the amplifier used to control the TEM. This amplifier is controlled in high-gain feedback, therefore it adjusts its output voltage to compensate for the V Sm that acts as a back EMF type voltage. This voltage V Sm in V is known as the Seebeck e↵ect, and it generates a voltage based on the temperature gradient over the TEM.
Current I Figure 4: Illustration of the voltage profile following a current step I ( ). Since the current amplifier is in closedloop, the voltage ( ) is increased to compensate for the back EMF voltage V Sm .
Time constants. Solving (14) for 2 unknowns, i.e., S m (T avg ), R m (T avg ), is generally not possible. However, as suggested in [24], V Rm and V Sm manifest in di↵erent time scales. This di↵erence in time constants is illustrated in Figure 4. At time t a , a current command I is applied to the amplifier, causing an instantaneous step in electric potential V Rm . While V Sm only manifests after a sufficient time has passed and a thermal equilibrium is reached at time t b yielding a T = T h T c over the TEM. By explicitly exploiting this di↵erence in time constants, both S m (T avg ) and R m (T avg ) can be determined from (14) using voltage measurements.

Experimental identification setup
In this section the temperature-dependent parameters are identified using the approach presented in Section 4.1. A dedicated TEM parameter identification setup is designed to isolate the TEM from external influences and facilitate accurate estimation of temperature-dependent parameters In Figure 5 a schematic representation of the setup is shown. The TEM is clamped between two stainless steel blocks to provide some additional thermal mass and spread the heat evenly. On the top, the hot side, the steel block is conditioned using a water cooling block and water chiller to provide a controllable temperature stable heat sink.  The setup is encapsulated by a 3D-printed enclosure made of High Impact PolyStyrene (HIPS) that is printed with a low infill of 10% to provide thermal insulation from the environment. The temperature measurements are obtained by using NTCs. Each stainless steel block contains 2 sensors, as indicated in Figure 5, where T 2 and T 3 are considered the TEM hot and cold side respectively. To mitigate heat transfer from the setup to the enclosure, small tabs connect the lower block to the HIPS enclosure, as shown in Figure 6, to minimize the contact area.

temperature-dependent identification
In this section, the method proposed in Section 4.1 is utilized to estimate the temperaturedependent parameters R m (T avg ), S m (T avg ) for multiple TEMs. To yield an accurate model for the purposes of this paper, a significant temperature range for T avg must be considered. To achieve this, the input current I is changed in small increments covering a wide operating range, as shown in Figure 7a.
Identification procedure. The identification procedure of the temperature-dependent parameters is described in Algorithm 1 where I are the steps in the current reference for the amplifiers, as shown in Figure 7a and I 0 = 0 is the initial current. The electrical resistance R m (T avg ) is estimated from the instantaneous voltage jump V Rm , shown in Figure 7b, that occurs after a step in current, since R m (T avg ) = V Rm /I. Then, the current is maintained until the system reaches a steady-state and associated T , as shown in 7c. This yields a back EMF voltage due to the Seebeck e↵ect V Sm = S m (T avg ) T , that is used to estimate S m (T avg ). By repeating this process both parameters are estimated for a range of T avg .
Temperature-dependent parameters. The identification procedure is repeated for 3 di↵erent, but of equal type, TEMs. This allows the characterization of an average parameter over a batch of actuators. While individual calibration curves could yield superior results, most applications do not allow for dedicated unit calibration tests, since they are both time intensive and expensive. The results of the identification procedure are shown in Figure 8 for R m (T avg ) and in Figure 9 for S m (T avg ). Both parameters show a linear dependency on T avg and a significant change of their value over the range of 15 [ C] to 55 [ C], that corresponds to a coldside temperature range of 5 [ C] to 85 [ C] since the hot-side is maintained at a constant temperature of 25 [ C]. The results show a small spread over the di↵erent TEMs, however, a larger sample size is required to yield a more definitive variance metric on the physical parameters. Moreover, since the input current profile consists of both positive and negative I steps, at similar T avg , some insight into possible hysteresis e↵ects is gained. In the results, the positive and negative current perturbation yield similar parameter estimates, indicating that hysteresis e↵ects are negligible.

Experimental validation
To verify the e↵ectiveness of a model structured as in (8) with temperature-dependent parameters     S m (T avg ), R m (T avg ) a comparison of modeling with and without temperature dependencies is shown in Figure 10, partially fulfilling Contribution C5. The modeling error, without including temperature dependencies as shown in Figure 10a, is quite significant, especially when a large temperature range is considered. By including the temperature dependencies of the model parameters the modeling error, as shown in Figure 10b, is significantly reduced. More importantly, while there is still a residual error, it appears to be much less correlated with the temperature range. This indicates that the temperature dependency is successfully taken into account.
G lin x Figure 11: Closed-loop control diagram including the feedback linearization law that linearizes the non-linear system G nl to obtain a system G lin that is linear in input v to output y = xc dynamics.

Feedback linearization
In this section an approach is presented that incorporates temperature-dependent parameters into a feedback linearization approach, with a stability guarantee. Thereby providing Contribution C3. The model obtained in Section 3 is leveraged to improve accuracy of the linearization over a wide operating range.

Feedback linearization
In contrast to earlier results in [5], the method in this paper considers temperature-dependent parameters in the linearization process. Consequently, the linearization is valid over a wider temperature range, enabling implementation for a larger class of systems.

Remark 3.
To facilitate the presentation, in this section a single TEM is considered. For multiple TEMs, e.g., 2 modules as shown in Figure 3, the cold side state vector is x c 2 R n M with n M the number of TEMs in the setup.
Given the system in (8) with the non-linear function F NL in (9). Consider the feedback linearization law I = (v, x), as shown in Figure 11, that provides I given a variable v and states that can be solved for I = (v, x) to find where S m (T avg ) and R m (T avg ) are both a function of the average temperature T avg = x h x c 2 .
Clearly, solving (15) for I has 2 possible solutions with (16) being the one that yields to smallest I for a given v. The feedback linearization law (15) is constructed for x c and not x h since the cold side of the TEM is most commonly the side of interest. A similar approach could be used to linearize towards x h , although the stability proof will be slightly more involved. The variable v has a clear physical interpretation, it represents the e↵ective heating power in the TEM. Moreover, a lower bound on v is required to ensure non-complex values in (16) given by Finally, an optimal current that maximizes the effective heating power is given by that is obtained if the e↵ective heating power (15) is minimized for I. Applying the feedback linearization law (v, x) to the system in (8), with the structure as proposed in (10), yields the following state equations Here, the "hot side" dynamics are non-linear since they are a↵ected by f NL h (x h , x c , u(v)) = 1 2 R m u 2 +S m x h u+K m (x c x h ), and the "cold side" dynamics are linear since they are solely a↵ected by v. By including the conduction term K m (x h x c ) in (v, x) a uni-directional coupling is achieved, i.e., the hot side does not influence the cold side. Consequently, classical linear control can be employed to generate an appropriate control input v for the linear plant G lin where the dynamics from v to y = x c are now linearized. The controller C lin is typically a collection of linear filters, e.g., a PID controller, that can be tuned through loopshaping techniques [25].

Stability
The feedback linearization law from Section 5.1 achieves unidirectional decoupling of the hot side to the cold side. Moreover, a controller for G lin can be constructed such that it is stable. However, this does not guarantee that an unreachable set point does not induce unstable behavior, known as thermal runaway. This can occur when an infeasible amount of cooling for v is required to achieve a certain temperature setpoint that is not achievable by the TEM. To guarantee full system stability, input to state stability can be guaranteed by constructing a Lyapunov candidate function where P H = P T H 0 and P C = P T C 0 are symmetric and positive definite matrices. The derivative of V in (20) is given, after explicitly exploiting the structure in (19), bẏ , assuming a solution exists that allows the structure for P H , P C as illustrated in (20). To guarantee the structure in (20) additional constraints are enforced during the solving of the LMIs in (23). Input to state stability can be concluded iḟ V < 0 in (22) and a solution to (23) exists under the structure imposed by (20). To conclude thatV < 0 holds, (22) is further simplified by employing the identities where min (Q) denotes the minimal eigenvalue of Q and x H can be separated into x L h , x h by using Figure 12: temperature-dependent bounds displayed as a surface for varying hot and cold side temperatures. Here, a positive current provides a flux from Tc to T h , and it can be observed that when Tc is high, more cooling is allowed.
to finally yieldV Since min (Q) > 0 for a positive definite matrix Q 0 and the ambient air w 0 is assumed to be bounded, the first termV 1 can be assumed to be strictly smaller than 0. To conclude stability by ensuring thatV < 0 it must hold thatV 2 < 0.

temperature-dependent stability bounds
The bound onV 2 < 0 in (31), is solved for I to yield the following temperature-dependent bounds on the allowable TEM current I [A] S m (T avg ), R m (T avg ) are the temperaturedependent parameters. The bounds on v [W] can be derived using a similar approach.
The bounds obtained in (32) are typically incorporated into the feedback linearization as a saturation bound on the TEM current. However, they can be included in a more advanced type of control, e.g., model predictive control, to apply the maximum allowable current given measured hot and cold side temperature. In Figure 12 the bounds in (32) are shown for varying T h , T c since they are temperature-dependent. Moreover, the accuracy of the bounds is improved by including temperaturedependent parameters S m (T avg ), R m (T avg ). In and a varying T c , both for constant parameters S m (308), R m (308) and temperature-dependent parameters S m (T avg ), R m (T avg ). The results show that the bounds vary for di↵erent values of T avg = T h +Tc 2 , and at T h = 318 [K], i.e., T avg = 308 [K], the bounds are equal.
Remark 4. The bounds in (32) can be conservative, as they form a su cient but not necessary condition for stability. The amount of conservatism depends on the specific solution to the LMIs in (23).

Experimental validation
By employing the high-fidelity model obtained in Section 4 combined with the procedure provided in Section 6 the application setup is linearized. Leveraging the feedback linearization law (v, x) in (16) the system G lin , as shown in Figure 11, is now linear To verify the e↵ectiveness of the feedback linearization, a comparison of the systems step response with and without linearization is made for di↵erent step magnitudes, contributing to Contribution C5.
Without linearization. The results without feedback linearization are shown in Figure 14, it shows both the step responses of di↵erent magnitudes and the normalized responses. The normalized responses in Figure 14b clearly show that the response varies with step magnitude. This indicates that the TEM system does not abide by the scaling principle of linear systems since the system is both state dependent and non-linear in the input.
With linearization. The results with feedback linearization, i.e., system G lin , are shown in Figure 15. Again, step responses of di↵erent magnitudes are shown in Figure 15a that are then normalized and shown in Figure 15b. The normalized responses, and especially y 2 , show significantly reduced correlation with the magnitude of the step. This illustrates that the system now behaves linearly in its input to output dynamics, thereby validating the e↵ectiveness of the feedback linearization. Moreover, by including the temperature dependency of S m , R m in the linearization law is now valid over a significantly larger temperature range when compared to earlier results [5].
Remark 5. Alternatively, it is possible to analyze the extent of the non-linear residual after feedback linearization through an analysis as described in   (b) Normalized step responses, and corresponding zoom plot of the initial response, using feedback linearization for both module 1 y1 and module 2 y2. Figure 15: Normalized experimental step response with feedback linearization. The normalized step response is almost completely invariant to the amplitude of the step, illustrating that the non-linear system is now linearized in input to output dynamics. [26]. This is done using multiple experimental realizations to determine the Best Linear Approximation (BLA) of the transfer function G lin in Figure 11 accompanied with the level non-linear distortions.

Observer design
In Section 5 an approach to feedback linearize the TEM based system is presented. This technique requires state information that is not always available. Moreover, in industrial applications it is often unwanted and sometimes impossible to attach numerous sensors to the TEM. Therefore, in this section, a temperature-dependent observer based approach is presented. It expands on existing results in literature [12] by explicitly accounting for temperature-dependent parameters. By combining a limited number of sensors with an internal model a state estimate is constructed, thereby providing Contribution C4.

State Estimator
Consider again the system in (8) but now the linear conduction terms K m (x c,h x h,c ) are combined into the A matrix. The mass matrix E can be taken to be unity or incorporated into the A and B matrices and f NL by multiplying those contributions with the inverse of E. Moreover, the system is now described in discrete time by an appropriate transformation, e.g., Euler approximation. This yields 2 6 6 4 as the state-space equations. Here, 2 R m u 2 and ↵ = T s a scaling factor that is equal to the sampling time T s for Euler approximation. The matrix C d relates the outputs to the states. Moreover, it is equal C d = C to the output matrix in continuous time and dependent on the sensor placement in the setup.
To construct a state-estimator an additional innovation term is added to (34) to yield, with slight abuse of notation, where L is the observer gain andx the state estimate. The observer gain is designed such that the estimation error,x(k) = x(k) x(k), dynamics are stable and converge within reasonable time. For a linear stable system, this can be always be achieved if the states can be reconstructed in the output of the systems, i.e., if they are observable. The observability of the states is confirmed, assuming the system is stable, using the Hautus test on the pair The estimation errorx dynamics are given by 2 where A o are the linear error dynamics and F NL o the non-linear contribution that includes the statedependent elements of f NL in (35).

Observer gain
The observer in (35) is constructed as a linearquadratic estimator (LQE), i.e., a Kalman filter. To determine the observer gain L that projects (y(k) C dx ) onto the state estimatex a discrete algebraic Ricatti equation is solved. Given a discrete system and noise covariance matrix of the process and measurement noise respectively given by the gain matrix L can be derived by solving a discrete-time algebraic Riccati equation assuming that v(k) and (w) are uncorrelated, i.e., E(w(k)v(k) T ) = 0.

Stability
In the typical application of the LQE, stability of the estimation errorx dynamics is ensured by evaluating the eigenvalues of A o = A d LC d . However, in the current application this is not su cient since F NL o also influences the estimation error as seen in (36), making the system non-linear.
To verify stability of the observer structure employed in this paper, a Lyapunov cadidate function is constructed as with P 2 R nx⇥nx 0, P T = P . The discrete time derivative of V (x(k)) = V (x(k + 1) V (x(k)) is given bỹ that must be negative for allx to ensure stability. This is equivalent to Clearly, the matrix F NL o is a function of S m (T avg ), which is temperature-dependent, and the input u.
To provide stability against a worst case F NL o a lower bound on F NL o is calculated by using the inequality Using this inequality it is possible to determine the worst case F NL o by determining the maximum eigenvalue max (F NL o ) over a range of T avg and u. Moreover, since F NL o is diagonal, see (36), its eigenvalues are determined by the entries on the diagonal itself. Therefore, the maximum eigenvalue is obtained using S M (T avg ) with T avg = 328 [K] Celsius, where S m is largest as shown in Figure 9, and u = 4 [A] which is well above the admissible current for the TEM. This then yields a worst case scenario bound for F NL o that is incorporated in (44). Stability of the estimation errorx(k) is then achieved if a solution is found to the set of LMIs in (44) using F NL o with the worst case parameters as determined by the lower bound in (46). If the LMIs in (44) can not be solved for a given observer gain L, stability can not be guaranteed and the design should be reconsidered by redesigning L. Remark 6. Analogous to the stability of the feedback linearization in Section 5, a solution to the LMIs (44) is only a su cient and not a necessary condition for stability of the estimation error dynamics.

Experimental validation
In this section an observer, as described in Section 6, is used to form a state reconstructionx using a limited set of sensors on the setup shown in Figure 1, completing Contribution C5. More precisely, only sensors on x 6 and x 7 , as shown in Figure 3, are used and the other sensors are used to evaluate the state reconstruction error.
In Figure 16 the results for 2 separate experiments are shown. For both module 1, with results in Figure 16a, and module 2, with results in Figure 16b, a sequence of step excitations is used. In the bottom half of the figures the reconstruction errorx = x x is shown, where x are the measured states andx the reconstruction made by the observer. It can be seen that the reconstruction error is quite small over the full temperature range. Moreover, the results are in the same order of magnitude as those shown in Figure 10b, indicating that the reconstruction by the observer is accurate up to the accuracy of the underlying model.

Conclusion
Advanced thermal control is a crucial area of research and development, especially in the medical, high-power lighting, and semiconductor industry. Controlling thermal dynamics using heater based actuators is inherently limited since they can only provide a positive heat flux. Using thermoelectric elements alleviates these limitations by providing active temperature control capable of both heating and cooling. This paper presents a full framework for temperature-dependent modeling and control of thermoelectric elements. Using temperaturedependent parameters yields a high-fidelity first principles model suitable for control. Using feedback linearization combined with an observer facilitates advanced control with a limited number of temperature sensors. It is shown that the framework yields accurate models and reliable control strategies by explicitly incorporating the temperature dependency of the model parameters, making it suitable for a large operating temperature range. This makes the framework readily applicable to a wide class of systems.  K m1 (x 11 x 9 ) + 1 2 R m1 u 2 + S m1 x 9 u K m2 (x 12 x 10 ) + 1 2 R m2 u 2 + S m2 x 10 u K m1 (x 9 x 11 ) + 1 2 R m1 u 2 S m1 x 11 u K m2 (x 10 x 12 ) + 1 2 R m2 u 2 S m2 x 12 u 0 1⇥6 Full matrices for the model given in (8) to facilitate reproduction of the results. Note that to facilitate the presentation, some significant digits are truncated from the matrices that might influence the results. In this paper, the full matrices are used during all simulations and experiments. Rob van Gils received the M.Sc. and Ph.D. degree ('08 and '12) in the Dynamics and Control group (Mechanical Engineering) from the Eindhoven University of Technology, Eindhoven, the Netherlands. For his PhD he focused on the feedback control and modelling of a thermo-dynamical system. Currently, he is a competence leader and senior technical specialist for the Thermal & Flow group with Philips Innovation Services, Eindhoven, The Netherlands.
Bram de Jager received the M.Sc. degree in mechanical engineering from Delft University of Technology, Delft, The Netherlands, and the Ph.D. degree from Eindhoven University of Technology, Eindhoven, The Netherlands. He was with Delft University of Technology and with Stork Boilers BV, Hengelo, The Netherlands. He is currently with the Eindhoven University of Technology. His research interests include robust control of (nonlinear) mechanical systems, integrated control and structural design, control of fluidic systems, control structure design, and applications of (nonlinear) optimal control. His research interests are in the field of data-driven modeling, learning, and control, with applications in precision mechatronics.