Soft Sensing-Based In Situ Control of Thermofluidic Processes in DoD Inkjet Printing

This article introduces a closed-loop control strategy for maintaining consistency of liquid temperature in commercial drop-on-demand (DoD) inkjet printing. No additional sensors or additional actuators are installed in the printhead while achieving consistency in liquid temperature. To this end, this article presents a novel in situ sensing-actuation policy at every individual liquid nozzle, where the jetting mechanism has three distinct roles. It is used for jetting liquid droplet onto the print media based on the print job. It is used as a soft sensor to estimate the real-time liquid temperature of the jetting nozzle. While not jetting liquid, it is used as a heating actuator to minimize the gradient of liquid temperature among nozzles. The soft sensing-based in situ controller is implemented in an experimentally validated digital twin that models the thermofluidic processes of the printhead. The digital twin is scalable and flexible to incorporate an arbitrary number of liquid nozzles, making the control strategy applicable for future designs of the printhead.

Soft Sensing-Based In Situ Control of automated fashion [2]. Since then, inkjet-based digital printing has been improved to be a faster, more cost-effective, and energy-efficient process. The noncontact and additive nature of inkjet printing has also attracted versatile applications. With the emergence of additive micromanufacturing, inkjet printing is now applicable to a wide range of materials, including polymers and metals, and also various media such as textile, wood, and circuit boards [3]. Due to the technical developments in inkjet-based digital printing, we are now gradually transcending "from the world of paper printing to printing the world" [4].
Today, the printer is a commercial peripheral, which produces printed products based on the user's demand. Due to the constant increase in the market competition of inkjet printing, the print quality stands out to be a decisive factor. In general, the print quality indicates how accurate the printed product resembles its digital counterpart. However, the quantitative metric that determines the quality of the printed product depends mainly on the application area and the user group. For example, in paper printing, the measure of quality is typically determined by the reflectance of the printed area, print density variations, offsets of the printed pixels, and so on [5]. For textile printing, the measure of quality is the gloss of the printed surface, texture geometry, surface smoothness of the printed area, and so on [3].
As the application domain of inkjet printing is quite diverse, printer manufacturers concentrate on how to improve the general process of inkjet printing instead of considering any application-specific metric. To this end, the primary motivation behind this research is to enhance the properties of the liquid droplets to improve the print quality, irrespective of the printing medium and its applications. One can view the inkjet printing process as a physical integration of liquid droplets with a solid medium. This process is also known as jetting. Based on the user-defined print job (a set of to-be-printed images), the jetting process deposits droplets of liquid on a solid medium (e.g., a sheet of paper). A set of printheads is responsible for generating and jetting the droplets in the right sequence, at the right position, with a prespecified volume and with well-defined (or carefully controlled) physiochemical properties.

A. Drop-on-Demand (DoD) Inkjet Printhead
A DoD inkjet printhead jets droplets of liquid based on the specific image that the user demands to print. A typical architecture of a commercial DoD printhead is shown in Fig. 1(a) (for instance, see XAAR printhead [6]). Here, the bottom stage of the printhead is called the nozzle platform (NP), which is divided into two mutually insulated parts NP left and NP right . Each part of NP consists of an array (divided into rows and columns) of nozzles. An example of the nozzles' placement with respect to an NP is shown in Fig. 1(b). The NP can be removed from the printhead and replaced with a new design. This feature allows for accommodating arbitrary numbers of nozzles without redesigning the entire printhead. On arrival of a print job, the following aspects play key roles in the jetting process.
1) DoD Bitmap: The user-defined image is translated into a bitmap. A bitmap stores the sequence in which nozzles are used for drop formation.
2) Liquid Inlet and Liquid Recirculation: Based on the image, the required flow rate for the deposition and recirculation of liquid is predetermined and assigned for every individual nozzle.
3) Mechanism for DoD Drop Formation: Every individual nozzle comprises a conduit, a liquid chamber, and a mechanism for ejecting liquid from the chamber, through the conduit, onto a printing media. This jetting mechanism differs based on the type of inkjet printer. For thermal inkjet printers, the jetting mechanism involves a resistive circuit, whereas another type of inkjet printers is equipped with a piezoelectric element attached to a diaphragm. In this article, the focus will be on the printers that are equipped with piezoelectric element. Depending on a specific bit map, a voltage pulse is applied to the piezoelectric element using a drive circuit. Actuation of the piezoelectric element forms a pressure wave inside the nozzle to create a droplet of liquid with predefined volume and push it out onto the medium. The frequency at which the jetting mechanism is actuated is in the range of kilohertz [5, pp. 9-10].

B. Quality Limiting Aspects in DoD Inkjet Printhead
In DoD inkjet printing, print quality depends on the physiochemical properties of every individual droplet (these properties are viscosity, surface tension, temperature, and the speed of sound inside the nozzle). During printing, these properties of droplets may deviate from their desired values. As a consequence, the predefined actuation pulse does not jet the droplet of desired physiochemical qualities and, therefore, results in poor print quality.
Two factors play crucial roles behind the deterioration of the droplet properties: 1) nozzles acoustics and 2) the temperature of liquid droplets. In [7], the problems related to nozzle acoustics are addressed using a feedforward controller. The current article considers that adequate measures are already incorporated to consider the problems related to nozzle acoustics and, hence, will not be discussed. On the other hand, during the jetting process, temperature gradient among the nozzles 1 may occur due to the following reasons.
1) Heat Dissipation in Nozzles: During jetting, piezoelectric element dissipates heat [5, p. 2]. Moreover, depending on drop volume, a jetting droplet results in impeding thermal energy. In particular, if not all the nozzles in the NP are simultaneously used for jetting, a gradient in liquid temperature among jetting and nonjetting nozzles takes place [8, p. 477].
2) Thermal Crosstalk in Solid and Liquid: In Fig. 1(b), the array of nozzles is physically coupled via the solid structure of the NP. Due to the different physical properties of solid and liquid, their mutual crosstalk causes changes in temperature among nozzles as well as adjacent printheads.
Due to such gradient in liquid temperature among nozzles, once a preset jetting pulse is applied, the resulting liquid droplet's temperature is not maintained to its desired value causing deterioration of the droplet's physiochemical properties. In general, the gradient in liquid temperature among nozzles occurs due to the exchange of thermal energy (heat flux) on mutually interacting components. The exchange of thermal energy (heat flux) on mutually interacting solid and liquid substances is a thermofluidic process. Concerning inkjet printing, the effect of this thermofluidic process is the gradient of liquid temperature among nozzles, and poor print quality is the immediate consequence.

C. Contribution
In the printing industry, to improve print quality, the majority of the works focuses on better understanding the process of drop formation or compensating nozzle acoustics, not specifically the effect of thermofluidic processes. For example, physical modeling and experimental methods of drop formation and relevant fluid dynamics in inkjet printing are studied in [2], [9] and [10]. The applications of control techniques are largely restricted to designing the DoD voltage pulses for generating droplets and compensating problems related to nozzle acoustics (cf. [5], [7], [11]). Apart from these works, there are only a few openly accessible patents that focus on thermofluidic processes, although they regard thermal inkjet printers and not piezoelectric ones. These methods conceptualize insulating apparatus and driving circuits that control the temperature of liquid ink by means of novel hardware design (see [12]). Hence, the state of the art for controlling thermofluidic processes in inkjet printing is still restricted to the hardware level. In contrast, a system-theoretic outlook to model and control the thermofluidic aspects is still mostly missing as far as DoD inkjet printers are concerned.
This article introduces, for the first time, a model-based feedback controller for the DoD inkjet printhead that uses no additional sensors and no additional actuators to compensate for the fluctuation of liquid temperature among nozzles. The applicability of the control system is demonstrated with the help of three novel contributing aspects.
1) The modeling framework presents a digital twin of the printhead, which is modular toward an arbitrary number of nozzles. Changes in the number of nozzles do not require redevelopment of the digital twin from the outset. 2) The piezoelectric elements that are already installed in every individual nozzle are used as collocated soft sensors. The self-sensing capability of a piezoelectric element is exploited to develop a data-driven algorithm and estimate the liquid temperature at every individual nozzle. 3) A concept of in situ actuation is introduced to control the fluctuation of liquid temperature. Here, a sensing-actuation policy is developed such that the controller uses the bitmap to anticipate the change in temperature at every nozzle and utilizes only the nonjetting nozzles to compensate for the temperature gradient.

D. Organization of This Article
The remainder of this article is organized as follows. After a brief explanation of the notations in Section II, a graph-theoretic modeling framework is presented in Section III to design the digital twin. In Section IV, the piezoelectric element is used as a soft sensor to estimate the liquid temperature at every individual nozzle. In Section V, the developed digital twin and the soft sensor are validated using an experimental setup. A model predictive controller for compensating the fluctuation of liquid temperature is developed and demonstrated in Section VI. Section VII provides some concluding remarks.

II. NOTATIONS
The set of natural numbers is denoted by N. Given two positive integers 0 < n < m, the subset N [n,m] := {n, n + 1, . . . , m} ⊂ N denotes all natural number between n and m, including n and m. The set of n × n positive semidefinite and positive definite matrices is denoted by S n 0 and S n 0 , respectively. Let a be a complex number defined as a := x + j y with j = √ −1. Then, Re(a) denotes the real part of a, i.e., xm and Im(a) denotes the imaginary part of a, i.e., j y. The complex conjugate of a is defined by a * = x − j y. The Hermitian transpose of an m × n complex-valued matrix A is denoted by A H . For a matrix A ∈ R m×n , A † denotes the pseudoinverse of A. The array [ a b ] is often written as col(a, b). Remark 1: Due to confidentiality reasons, in this article, all the plots are normalized by subtracting the data points with a fixed nominal value.

III. DIGITAL TWIN OF THERMOFLUIDIC PROCESSES:
PHYSICAL MODELING The first objective of this article is to develop a digital twin of the jetting process for describing its thermofluidic aspects.
The digital twin constitutes a generic physical model of the thermofluidic processes that can be used throughout the design cycle. The digital twin must meet the following properties.
1) Modularity: The model is modular with respect to an arbitrary number of nozzles. 2) Flexibility: The model is flexible such that, during any stage of prototyping, necessary adaptation and modification on the model can be easily incorporated. 3) Versatility: The model serves as a versatile tool that can be used for simulation, design optimization, fault diagnosis, control design, and experimental validation. To cope with the ever-increasing demands on throughput, the modularity of the digital twin offers the freedom of increasing the number of nozzles without the need of remodeling. A flexible modeling approach can easily incorporate necessary adaptation for increasing accuracy or providing additional functionality.

A. Abstract Definition of Graph-Theoretic Thermofluidic Processes
This article presents a graph-theoretic framework that makes the digital twin modular, flexible, and versatile. The graph-theoretic model consists of the following three definitions.
1) Topology: A finite and connected graph consists of a set of nodes N and edges E. It is defined by the following triad: a) The time axis is defined by set T. If the thermofluidic processes are modeled in continuous time, T := [0, ∞). If discrete time is considered to model the thermofluidic processes, T := {kt d | k ∈ N}, where the sampling period t d is a fixed scalar. The time index is denoted by t ∈ T. b) Each object N i in the set N , i ∈ N [1,m] , is a node or component of the graph. c) The set of edges E = {E i, j | for all (i, j ) with A i, j = 1} describes the interconnection between a particular node and its neighborhood nodes. Here, A ∈ N m×m is an adjacency matrix that has either zero or one as its elements depending on whether the nodes N i and N j are physically attached or not. Precisely, its entries are 2) Edges: Every individual edge E i, j ∈ E is defined by the following couple: a) Wherever N i is interconnected to N j , i.e., A i, j = 1, there are time-varying interconnection signals associated with the edge. Let these signals be b) The relation among interconnection signals l i, j and l j,i defines a subspace M i, j ⊂ L I i, j . 3) Nodes: Every individual node N i is associated with the following triple: a) There are three categories of time-varying signals that belongs to a space S i associated with every ii) The to-be-manipulated control signals c i : T → C i . iii) Associated with every individual node, the collection of all interconnection signals is l i : T → L i . They are related to an individual edge such that L I i, j ⊆ L i × L j . iv) The user-defined print job introduces a set of signals that affect the node behavior. They are given by Given a specific bitmap, the signals d i are restricted by a subspace B i ⊂ D i . c) The behavior of N i is captured by a relation among all the signals (x i , c i , l i , d i ) and defines a subspace P i ⊂ S i . Depending on the application, all the subspaces in this definition are assigned by the user. For instance, it can be either a real vector space or space of squareintegrable/square-summable signals. Section III-B specifies these subspaces for the application of inkjet printhead.

B. Specification of Graph-Theoretic Thermofluidic Processes for a Particular Printhead
To implement a digital twin of the DoD inkjet printhead, the user simply requires to provide concrete specifications for each of the definitions (1)-(3). These specifications are given in the following items.
1) Printhead Topology: The digital twin of the printhead is a finite and connected graph according to the definition (1). 1) Each node represents a solid or liquid component. For a specific design, the user specifies the number of nodes. 2) The adjacency matrix A is defined based on the architecture of the printhead. 3) Except for the adjustable NPs, the remaining architecture of the printhead is typically kept identical. Then, specifying the arrays of nozzles is sufficient for defining the entire topology of a newly designed printhead. 4) The user defines T to specify whether the signals are in continuous time or discrete time. For the printhead shown in Fig. 1(a), individual nodes (components) and the interconnection topology are shown in Fig. 2.
2) Interconnection Among Components: The interconnection captures: 1) the communication; 2) sharing of information; and 3) exchange of (thermofluidic) energy among is an edge describing convective thermal energy due to the inflow of liquid along the inlet channels.
is an edge describing convective thermal energy due to the recirculation of liquid along the return channels.
interacting components. Concerning interconnected systems, the role of interaction is widely studied in [13]- [16]. In case of thermofluidic processes, the exchange of information and energy is either conductive (between two solid components) or convective (between a solid and liquid component or between two liquid components). Typically, a bidirectional exchange occurs between two solid or a solid and a liquid component. A unidirectional exchange occurs between two adjacent liquid channels toward the direction of the flow. Hence, a distinction is made between an edge E i, j and another edge E j,i to distinguish between bidirectional and unidirectional interconnections.
To model an individual edge between adjacent components N i and N j , i.e., wherever A i, j = 1, the user must specify the following items according to (2).
1) In terms of temperature of the adjacent nodes ( • C), the interconnection signals between N i and N j are cat- Here, the user specifies the dimension of the inter- 2) Based on the conductive or convective thermal energy exchange, point-wise in time, the subspace M i, j defines either unidirectional or bidirectional interconnection relations according to Wherever A i, j = 1, the user specifies the relation (4) by defining a constant matrix M i, j of suitable dimension.

3) Solid or Liquid Components:
For individual nodes, based on the physical properties and dimensions, the Biot number [17] justifies (experimentally assessed to be less than 0.1, e.g., individual nozzle's Biot number is found to be below 0.09) the possibility to consider a lumped model (finite dimensional) while neglecting the spatial variation of thermofluidic aspects in the component. A lumped model describes the evolution of temperature at individual node due to: 1) the conduction between two solid components; 2) convection between two liquid components; and 3) convection between a solid and a liquid component.
To model an individual node N i , the user must specify the following items according to (3).
1) The user specifies the dimension of the internal state variable x i (t) ∈ R n i x of individual node N i . In this case, x i is the temperature ( • C) of the solid or liquid component.
2) Control signals are categorized into control inputs (in terms of thermal power (watt) that can be manipulated) u i (t) ∈ R n i u and measured outputs (in terms of temperature ( • C) that can be measured) y i (t) ∈ R n i y . They define the control signals c i := col(u i , y i ).

3) Grouping together all the edges corresponding to an individual node, interconnection inputs
Based on the user-defined relation (4), these signals are constructed according to Using an user-defined bitmap, the flow parameters of inlet liquid and liquid return (m 3 /sec) are known and stored in a time-varying matrix (t) for all liquid channels. To allocate flow parameters for every individual node, an operator is the allocated flow parameter for the node N i . The function i (t) relates print-job associated inputs (in terms of convective heat flux, watt) p i (t) ∈ R n i p and outputs q i (t) ∈ R n i q (in terms of thermal energy over unit volume of liquid, kg • C/m 3 ) according to the following algebraic relation point-wise in time: Here, (5) describes convective thermal energy transfer due to the temperature difference in the liquid between the inlet and the return sections that depend on the liquid flow rate. Using user-defined (t) and T i in (5), the signals d i (t) and the matrix i (t) are automatically defined for individual node. 5) The dynamic relations among defined signals Equation (6) is derived using the laws of mass and energy balance (see [18]). Based on the user's choice of T, (6) can be used in either continuous time or discrete time. Accordingly, on functions f : T → R n , the operator Q t is defined as

C. Equivalent Classes of Thermofluidic Model
In the previous exposition, the graph-theoretic representation (1)-(3) presents a paradigm for building and upscaling the thermofluidic model for a DoD printhead with arbitrary number of nozzles. Subsequently, a compact multi-input and multioutput (MIMO) representation may appear to be useful for design optimization or synthesizing an observer-based controller. By rearranging signals and performing algebraic operations, there are three ways to equivalently represent the graph-theoretic model. 1) By Stacking Node Signals: 2) By Eliminating Interconnection Signals: 3) By Eliminating Print-Job Signals: Among them, the behavior of the graph theoretic model (1)-(3) is equivalently described byP I in (7). The equivalence holds as (7) is a result of signal rearrangement in (1)-(3). In particular, (x, c, l, d) are obtained by stacking the signals [1,m] . In (7), the matrix M is constructed by combining the relation v j,i = M i, j w i, j whenever A i, j = 1 for all i, j ∈ N [1,m] to form one algebraic relation among v and w. Similarly, the matrix (t) = diag( i (t)) i∈N [1,m] . The functions f x , g w , g q , and g y are obtained by stacking f i x , g i w , g i q , and g i y , respectively, in a column over all i ∈ N [1,m] . Furthermore, for each m, n ∈ {x, v, w, p, q, u, y}, A mn , B mn , C mn , and D mn are obtained by diagonally stacking A i mn , B i mn , C i mn , and D i mn over all i ∈ N [1,m] . One can also view P II and P III as projections of P I onto the subspaces spanned by the signals (x, c, d) and (x, c), respectively. Representing such projections using (8) and (9) requires elimination of interconnection signals l(t) and print-job signals d(t). The next theorem provides conditions under which such elimination of signals is possible.
Theorem 1 (Model Equivalence Under Signal Elimination): Proof: Deriving (10) requires algebraic manipulations to uniquely express v in terms of (x, u, p) by using v(t) =Mw(t). The uniqueness of this expression requires the invertibility of (I − D wvM ). The same argument holds for deriving (11), where d is uniquely expressed in terms of (x, u). In this case, one requires the invertibility of (I − D 11¯ (t)) for all t ∈ T.

IV. DEVELOPMENT AND CALIBRATION OF SOFT SENSOR
For control, monitoring, and fault diagnosis of the thermofluidic processes, it is necessary to acquire real-time information about the liquid temperature. However, this must be achieved without incorporating any additional sensors. To this end, the piezoelectric elements that are located at every individual nozzle become useful. Apart from jetting liquid droplets, the self-sensing capability of the piezoelectric actuator (see [19], [20]) is calibrated to make a sensing device from which liquid temperature can be estimated at every individual nozzle.

A. Mechanism of Acoustic Sensing
A piezoelectric element works in two operating modes. For jetting droplets of liquid, the piezoelectric element is in actuation mode where a sequence of trapezoidal voltage pulses allows droplet formation and ejection [21]. Using the same piezoelectric element, in sensing mode, one can capture the change in pressure oscillation and the nozzle acoustics during and after the actuation [22]. The mechanism is also known as acoustic sensing and the signal measured is called acoustic signal (volt). The acoustic signal provides useful information about the dynamics inside a nozzle (see [2], [5]).
By exploiting this self-sensing-actuation capability, in [2] and [19], two different approaches are presented to simultaneously operate the piezoelectric element in both actuation and sensing mode. However, these approaches require power-consuming and expensive electronic hardware that makes it difficult for online implementation at every individual nozzle [5,Ch. 3]. Moreover, these measurements are typically used for fault diagnosis and monitoring of the nozzle acoustics. For cheaper signal processing and its online application as a temperature sensor, in this article, a time delay (t delay ) is introduced between the actuation mode and sensing mode of the piezoelectric element. Once the actuation pulse is applied, after the time delay, an external trigger switches its operation to sensing mode, and the residual part of the acoustic signal is captured. Fig. 3 shows this sequential actuation-sensing scheme in an individual nozzle.
To use the piezoelectric element at every individual nozzle as a soft sensor for estimating liquid temperature, this article uses the residual acoustic signal. As specific trapezoidal pulses can be designed on the basis of the desired droplet temperature [5,Ch. 3], an experiment is performed to jet liquid at temperature T 1 • C and T 2 • C with T 1 < T 2 . The corresponding acoustic sensing signals are measured and shown in Fig. 4. The difference in the measured signals indicates that the acoustic sensing signals can be used for estimating liquid temperature.
Let the finite samples of acoustic sensing signal be given as the set {s n | n ∈ N [1,N ] }, with the length of the dataset as 0 < N < ∞, and the sampling period is fixed as T s . The aim is to estimate the parameters α, ζ , ω, γ , and φ such that (12) optimally approximates the measured signal {s n | n ∈ N [1,N ] }. To this end, the following algorithm is developed.
Algorithm IV.1 Reconstruction of Acoustic Sensing Signal: 1) Model Hypothesis: Consider the following class of stable, second-order, single-input-single-output systems: z(n + 1) = A s z(n) + B s w(n), y(n) = C s z(n). (13) Here, for n ∈ N, z(n) ∈ R 2 , and w(n) ∈ R is the applied impulsive input. The output signal y(n) ∈ R is modeled as (12). 2) Estimation Algorithm: It is given in Algorithm 1. The proof of the correctness of the algorithm as well as the assessment of its accuracy is included in the supplementary material. Here, the reconstruction of (12) acts as a digital filter to improve the noise level on the measured signal. Furthermore, the estimation of α, ζ , ω, γ , and φ also offers an online algorithm to monitor the dynamic behavior inside every individual nozzle.
2) Estimation of Liquid Temperature Using Acoustic Energy: A next step would be to determine how individual parameter α, ζ , ω, γ , or φ relates to the change in liquid temperature in a particular nozzle. A series of experiments is performed to correlate each of these key parameters with liquid temperature. It is found that the liquid temperature has the closest relation with the acoustic energy that is defined as the squared 2-norm of the estimated y(n), y(n) 2 2 (in Volt 2 ; equivalent to the notion of energy assuming that the equivalent resistance is constant). One expects lower acoustic energy in a lower temperature and higher acoustic energy in a higher liquid temperature. Algorithm 1 Estimating α, ζ , ω, γ , and φ.

C. Calibration of Energy-Temperature Characteristics
To estimate liquid temperature, a characteristic relation between acoustic energy and liquid temperature needs to be established for an individual nozzle. This is achieved by performing the following three experiments consecutively.
1) Experiment 1: In the first experiment, three operating points are chosen for the liquid temperature. At a fixed operating point of temperature, the liquid droplets are jetted using the actuation mode of the piezoelectric element. After that, it is switched to sensing mode 20 times consecutively, and each time the acoustic sensing signal is measured (see Fig. 8). Every individual signal contains 100 samples of measured data. Using the measured data, the acoustic sensing signal is reconstructed using Algorithm IV.1, and its energy is computed. Based on this experiment, the energy-temperature curve is shown in Fig. 5. Despite the expected trend in  change of acoustic energy with temperature, the variation of energy estimate is significantly large. Therefore, estimating liquid temperature using an energy-temperature curve is still unreliable.
2) Experiment 2: The evaporation of liquid changes its viscosity and temperature [18] and therefore also its acoustic energy. In order to present the consequence of evaporation, an experiment, similar to Experiment 1 (see Fig. 8), is performed where four operating points of temperature are chosen. Liquid droplets are jetted at these operating temperatures by using 2000 consecutive jetting pulses. The last jetting pulse is, thereafter, followed by switching the piezoelectric element in sensing mode and measuring a sequence of acoustic sensing signals. Here, the aim is to investigate the evolution of acoustic energy over time. In Fig. 6, the estimated acoustic energy of the acoustic signal is shown over 1 s. Over time, the effect of evaporation results in an initial decay of energy and the decay rate is higher if the temperature of the liquid droplet is higher. This experiment implies that, due to evaporation, estimating the relation between temperature and acoustic energy depends on the timed sequence (i.e., t delay ) between jetting and sensing. In other words, the larger the difference between the time of jetting pulses and the time of measurements, the lower the acoustic energy. As a result, the estimation method based on this experiment has a poor signal-to-noise ratio.
3) Experiment 3: To circumvent the effect of evaporation from the energy-temperature relation, only those acoustic sensing signals should be considered, which are measured immediately after the jetting pulse. To this end, another experiment is performed where the liquid temperature is varied over four operating points. In contrast to the previous experiments where a series of acoustic sig-    nals is measured one after another, this time, jetting of the liquid droplet is directly followed by the sensing of the acoustic signal. This sequence of consecutive jetting and sensing is repeated 20 times (see Fig. 9). Using the measured data, the acoustic energy is estimated. Fig. 7 shows that the adapted approach of consecutive jetting sensing significantly improves the variance in the estimate of energy-temperature relation.
In Figs. 8 and 9, the designs of Experiments 1 and 3 are compared. The key difference between them is the sequence in which jetting actuation and acoustic sensing are performed. Based on the sequential jetting-sensing mechanism devised in Experiment 3, the energy-temperature characteristic curves are estimated and calibrated for every individual nozzle. In Fig. 10, energy-temperature curves are shown for four different nozzles in the printhead.

D. Estimation of Liquid Temperature Using Energy-Temperature Curve
Based on the energy-temperature curves, one can determine a parametric linear model relating the acoustic energy and the liquid temperature. For an individual nozzle, let the energy-temperature characteristic curve be modeled as Here, for i th nozzle, x i is the liquid temperature and φ i is the acoustic energy. The unknown parameters m i , c i ∈ R are obtained by fitting the respective energy-temperature curve. Every time a liquid droplet is jetted from a nozzle, thereafter, it can be followed by measuring the acoustic signal. The acoustic signal is modeled as (12) and its energy φ i is determined using Algorithm IV.1. Subsequently, the corresponding liquid temperature x i is obtained from (14).

A. Experimental Setup
To validate the digital twin as well as the soft sensor, an experimental setup is built, as shown in Fig. 11. The setup consists of two liquid vessels that are connected with a printhead that has 160 nozzles. By keeping two vessels at a fixed level, a constant recirculation flow through the printhead is established. To raise the temperature of the liquid inlet in the printhead, the upper vessel is equipped with a heater. For each side of the NP, two locations are selected to place thermocouples (in Fig. 12, they are called T l,1 and T l,2 for NP left and T r,1 and T r,2 for NP right ). These locations are kept in the close vicinity of the liquid nozzles to receive real-time information of its temperature. Based on the architecture of the printhead in Fig. 2, T l,1 and T l,2 are close to the liquid node nM and n1 on the NP left , respectively. Similar arrangement goes for T r,1 and T r,2 on NP right . Moreover, these thermocouples are used solely for validation purpose and are not allowed in the final product. Furthermore, at each side of NP, one nozzle is selected whose temperature is monitored using the developed soft sensor (in Fig. 12, they are called N l corresponding to the node nM for NP left and N r corresponding to the node nM for NP right ).

B. Experiment Design
The printhead is operated to print 3500 A-4 sheets of paper. Only the nozzles in NP left are actuated for jetting while keeping the nozzles in NP right idle. The bitmap is designed such that the following holds.
2) In the next 100 s, another 500 pages are entirely printed by using all the nozzles in NP left . 3) In the next 100 s, another 500 pages are left unprinted. 4) In the next 100 s, another 500 pages are entirely printed by using half of the nozzles that are located on the left side of NP left . 5) In the next 100 s, 500 pages are left unprinted. 6) In the next 100 s, another 500 pages are entirely printed by using half of the nozzles that are located on the right side of NP left .

7)
In the next 100 s, the last 500 pages are left unprinted.

C. Setting Up the Digital Twin for Simulation
Based on the bitmap, the corresponding flow parameters are assigned with the print-job signals. With 160 nozzles, a graph theoretical digital twin of the printhead is built by following the nodal structure that is shown in Fig. 2. Every individual node is associated with its temperature (in • C) as its internal state. Based on the topology of the printhead, the interconnection signals and their relations are built. The print job also assigns the flow parameters to every individual node through the print-job signals and their relations. Once the model is built using the definitions (1)-(3), anyone of the representations P I , P III , and P III can be used for simulating the digital twin.

D. Results
In Fig. 13, at location T l,1 and T l,2 and T r,1 and T r,2 , temperature data measured by the thermocouples are compared against the temperatures of the corresponding nodes simulated by the digital twin (according to Fig. 2, nM and n1 on NP left as well as NP right , respectively). Using the developed soft sensor, in Fig. 14, the estimated liquid temperature at N l and N r is compared against the temperatures of the corresponding nodes simulated by the digital twin (according to Fig. 2, nM on NP left and NP right , respectively). The developed model captures the thermofluidic behavior of the printhead by indicating the offset in temperature among jetting and nonjetting nozzles. For example, as the nozzles in NP right are idle during the entire print job, its temperature variation is significantly smaller. Nevertheless, there is still a small temperature variation for nozzles in NP right due to thermofluidic crosstalk among nozzles. Using only the left half of the nozzles or the right-half of the nozzles in NP left , among nozzles, there is a significant Fig. 13. Comparisons among the temperature measurements from the thermocouple and the digital twin. The legends are as follows: temperature at T l,1 ; temperature at the corresponding node in the digital twin; temperature at T l,2 ; temperature at the corresponding node in the digital twin; temperature at T r,1 ; temperature at the corresponding node in the digital twin; temperature at T r,2 ; temperature at the corresponding node in the digital twin. estimated temperature of N l by soft sensor; temperature at the corresponding node in the digital twin; estimated temperature of N r by soft sensor; temperature at the corresponding node in the digital twin. temperature difference over time (around 2 • C difference between T l,1 and T l,2 during 300-400 s).
Similarly, in Fig. 14, the soft sensor's estimation of the liquid temperature follows the same trend that of the digital twin. The sudden oscillations that appear in the soft sensor's estimation are due to sporadic bubble entrapment or drying of liquid inside the nozzle during jetting, which may fluctuate the liquid temperature. Such erratic phenomena are not captured in the digital twin. On the other hand, the nonjetting nozzles' temperature estimation is almost identical to the model.

VI. In Situ CONTROLLER: A PROOF OF PRINCIPLE
The performance limiting aspects of thermofluidic processes is the fluctuation in liquid temperature among adjacent nozzles in a printhead and among adjacent printheads. To maintain temperature consistency, a controller is synthesized that does not require incorporating any additional sensor or actuator. The only resources the controller uses are: 1) the bitmap as prior knowledge about the flow pattern and jetting sequence of nozzles in the printhead; 2) the model derived in Section III to anticipate and predict the evolution of thermofluidic behavior in the printhead; and 3) the soft sensor as feedback information on temperature and viscosity at every nozzle as derived and discussed in Section IV.

A. Concept of In Situ Sensing-Actuation
Every individual nozzle is equipped with a piezoelectric element for jetting droplets of liquid. By applying a voltage, its resistive property allows the piezoelectric element to dissipate heat. In this controller design problem, the resistive property of a jetting actuator is used as a source of heating. Thus, the piezoelectric element makes every individual nozzle equipped with a local, independent, in situ self-sensing heating actuator.
However, the piezoelectric elements can only be used as heaters for control input when there is no need for jetting liquid. To implement the in situ sensing-actuation scheme, at all time t ∈ T, every individual nozzle is in either jetting, heating, sensing, or idle mode. To allocate the nozzles in these four modes of operations, the following time-varying set-valued maps are defined.
1) Jetting nozzles: At time t, let there be n p (t) number of nozzles that are used for jetting. Hence, 2) Heating nozzles: At time t, let there be n h (t) number of nozzles that are nonjetting and used as heaters. Hence, 3) Sensing nozzles: At time t, let there be n s (t) number of nozzles whose temperatures are sensed. Hence,

4)
Idle nozzles: At time t, let there be n u (t) number of nozzles that are in none of the above modes. Hence, The cardinality of each set (as a function of time) depends on a specific bitmap. Intersections of S p (t) and S h (t) must be empty for all time t ∈ T. This requirement on functional exclusion is a serious engineering challenge as a piezoelectric element in heating mode is not supposed to jet. As described in Section IV, there is always a delay between applying the jetting pulse and sensing the temperature using a soft sensor. Therefore, the jetting nozzles are chosen not to be used for sensing (i.e., the intersection between S p (t) and S s (t) is chosen to be empty).

B. Performance Specifications
The controller must satisfy the following performance specifications 1) Irrespective of the bitmap, the differences in temperature among adjacent nozzles must be below ±0.3 • C.
2) The usage of thermal actuation must be limited to sustain the operational lifetime of piezoelectric elements. 3) Cooling of the liquid is not possible using thermal actuation. 4) Thermal actuation of piezoelectric material should be nonjetting, i.e., it should not form new droplets of liquid. Intuitively, using the in situ sensing-actuation scheme, the controller is expected to actuate only on a few heating nozzles [from the set S h (t)] adjacent to the jetting ones to limit the temperature difference among nozzles while satisfying the performance specifications. At the same time, the nozzles' temperature can be monitored using the implemented soft sensor at the sensing nozzles [from the set S s (t)].

C. Generating Voltage Pulse for Heating Actuators
To apply the required amount of thermal power as control inputs, one needs to design voltage pulses that are to be applied on the piezoelectric elements in heating mode. Such voltage pulse must not cause ejection of droplets while satisfying the power requirement. The principle behind generating a nonjetting and heating voltage pulse is to modulate the amplitude and frequency of the trapezoidal pulses (see [23]). Recently, in [24], a mechanism is developed so as to generate specific thermal energy while not interfering with the bandwidth that defines the jetting mode. Using bandpass modulation, the width and the height of the trapezoidal voltage pulses are tuned. Further details on these signal implementations are omitted for brevity.

D. Model-Based Prediction of Liquid Temperature in Nozzles
As mentioned in Section III, there are four equivalent ways to represent the model of thermofluidic processes and one has the freedom to use anyone of them for predicting the future evolution of liquid temperature in nozzles. Here, the thermofluidic behavior is considered to be represented according to P III in (9). In discrete time, i.e., for t ∈ {kt d |k ∈ N}, the behavior of P III is governed by the following equations: Here, Q t (x)(t) := x(t + t d ), t ∈ {kt d |k ∈ N}. If there are n nozzle number of nozzles in the printhead, u(t), y(t) ∈ R n nozzle .
During jetting, the actuation of jetting piezoelectric actuators causes heat dissipation [23]. Given the amplitude and frequency of the voltage pulses, one can experimentally estimate the thermal power (in Watt) generated in jetting nozzles over time, and they are modeled as known disturbances d(t) ∈ R n nozzle .
However, to implement the in situ sensing-actuation scheme, the bitmap allocates the nozzles according to the four modes of operation defined in Section VI-A. To map the total number of nozzles to its allocated modes of operation, binary matrices S p (t) ∈ R n nozzle ×n p (t) , S h (t) ∈ R n nozzle ×n h (t) , and S s (t) ∈ R n s (t)×n nozzle are introduced. Their (i, j )th entries are As a result, (15) is modified as

t)B(t)S h (t) S s (t)C(t) S s (t)D(t)S h (t)
. (16) Here, at time t ∈ T, the control inputs u h (t) ∈ R n h (t) are applied by the set of heating nozzles (in terms of thermal power, watt). The disturbance due to jetting d p (t) ∈ R n p (t) is applied by the set of jetting nozzles (in terms of thermal power, watt). The sensed outputs y s (t) ∈ R n s (t) are sensed liquid temperature by the set of sensing nozzles (in terms of temperature, • C).

E. Formulating Operational Constraints
Let u max ∈ R n nozzle be the maximum admissible thermal power of all piezoelectric elements. Here, u max is determined such that it takes the degradation of piezoelectric elements into account while preventing the ejection of droplets of liquid. On the other end, cooling is not allowed. As a result, at every time instant t ∈ T, heating piezoelectric elements must satisfy the constraint Having bounded actuation and bounded thermal disturbance during jetting, no additional constraints on the state variable x(t) have been enforced.

F. Specifying Control Criterion
The purpose of the controller is to control the temperature difference among nozzles. For the i th and j th nozzle that are adjacent to each other, the temperature difference is defined as z i j = x i − x j . Concatenating every individual difference of two adjacent nozzle's temperatures (i.e., column-wise stacking every z i j for each j ∈ { j |A i, j = 1} and every i ∈ N [1,n x ] ), one obtains the to-be-controlled variables as z(t) = H x(t), where z denotes the vector of individual temperature differences among two adjacent nozzles and H is a constant matrix.
At every time instant t = kt d with k ∈ N, the requirement is to control the liquid temperature gradient among nozzles over a finite horizon of future time instants t ∈ T k N , where To this end, a reference tracking problem is formulated to minimize the following cost functional: Here, the reference trajectories (x r , u r (t)) are predetermined from a steady-state model that equates the temperature difference of nozzles to zero. In other words, for t ∈ T k N , (x r , u r (t)) is the solution of the following linear equations: Moreover, for all t ∈ T k N , one must predefine the time-varying weights Q k (t) ∈ R n x ×n x and R k (t) ∈ R n h (t)×n h (t) in (18) to penalize the deviation of states and inputs from their respective references.

G. Formulating Model Predictive Control Scheme
A digital controller is required that minimizes (18) by using future prediction on the states in (16) while satisfying the constraints in (17). To this end, a reference-tracking model predictive control (MPC) scheme is presented.
At every time instant t = kt d with k ∈ N, for convenience, let any function f (t) be defined as f 0|k := f (kt d ). From its initial value f 0|k , over a finite horizon of future time instants . Starting from the time step kt d , the N-horizon MPC scheme amounts to minimizing the following cost functional (see [25]) (20) subject to the following constraints: Here, for i ∈ N [0,N −1] , For all t ∈ T k N , the reference trajectories (x r , u r ) are determined by solving (19). The weights Q ∈ R n x ×n x and R i|k ∈ R n h (t)×n h (t) are user-defined matrices. The terminal penalty Q k ∈ R n x ×n x in the cost (20) and F N |k , b N |k in the terminal constraint (21c) are specifically chosen to ensure the stability of the closed-loop system [26].
The above constrained optimization problem is solved for every k over the horizon kt d , . . . (k + N)t d of N future time samples. The sequence of predicted states and future inputs is denoted asx k := {x 1|k , . . . , x N |k } andū h k := {u h 0|k , . . . , u h N −1|k }, respectively. If the minimizer of the optimal control problem (20) and (21) is denoted byū h * k := {u h * 0|k , . . . , u h * N −1|k }, then its first entry u h * 0|k is implemented as control input at time step t = kt d . Subsequently, on a receding horizon, the optimal control problem is solved again at step (k +1)t d using x 0|k+1 = x((k + 1)t d ) as its initial state. In a receding horizon fashion, this procedure is continued iteratively over all time t ∈ {kt d | k ∈ N} [27].
One may substitute the equality constraint (21a) in (20) to eliminate the variables {x 0|k , . . . , x N |k }. This results in the following dense linearly constrained quadratic programming (LCQP) problem.
Optimization Problem 1 (Dense LCQP): Here, the decision variables are the to-be-applied control input U k = col(u h 0|k , . . . , u h N −1|k ). Moreover, X r k = col(x r , . . . , x r ) and U r k = col(u r 0|k , . . . , u r N −1|k ) are the reference values for states and inputs, respectively. Moreover, The following result summarizes the convergence properties of the resulting controlled system.
Theorem 2 (Stability of the In Situ Controller): Let the matrices R i|k ∈ S n h (t) 0 and Q ∈ S n x 0 are given for all time t ∈ T k N . Moreover, X k ∈ S n x 0 and Y k ∈ R n h (t)×n x satisfy the following linear matrix inequality (LMI): Furthermore, let (22) have a unique minimizer U * k := {u h * 0|k , . . . , u h * N −1|k } under the following conditions. 1) The terminal weight matrix Q k admits Q k = X −1 k . 2) The terminal constraint F N |k x N |k ≤ b N |k admits Fig. 15. Liquid temperature of 48 nozzles (equally distributed over two rows that correspond to the two rows at each subfigure) for Scenario 1 in controlled and uncontrolled case. denotes nozzles that are jetting. denotes the adjacent nozzles that are used as thermal actuators. denotes the nozzles that are not used for heating. (a) Controlled Scenario 1 (at 100 s): NP left is fully used for jetting. The maximum absolute difference of temperature between two adjacent nozzles is 0.1988 • C. (b) Uncontrolled Scenario 1 (at 100 s): NP left is fully used for jetting. The maximum absolute difference of temperature between two adjacent nozzles is 2.2551 • C.  denotes nozzles that are jetting. denotes the adjacent nozzles that are used as thermal actuators. denotes the nozzles that are not used for heating. (a) Controlled Scenario 3 (at 300 s): right half of NP right is used for jetting. The maximum absolute difference of temperature between two adjacent nozzles is 0.2 • C. (b) Uncontrolled Scenario 3 (at 300 s): right half of NP right is used for jetting. The maximum absolute difference of temperature between two adjacent nozzles is 0.9368 • C.
Then, setting u h (kt d ) = u h * 0|k , there exists δ > 0 for which every initial condition x(0) − x r 2 < δ and every solution to (16), k → x(kt d ), admits lim k→∞ x(kt d ) − x r 2 = 0. Proof: The proof is included in the supplementary material.

H. Implementation and Illustration of MPC Scheme
As proof of principle, the presented MPC scheme is implemented to minimize the difference in liquid temperature among nozzles. To this end, the following specifications are considered.
1) Printhead Configuration: The printhead chosen for implementing MPC has the architecture shown in Fig. 1(a). There are 48 nozzles that are equally divided over two sides of the NP, and at each side, 24 nozzles are equally distributed over two rows.

2) Bitmap Specification:
For generating a test bitmap, the following scenarios are considered for 300 s in total. 3) Setting Up the Digital Twin of the Model: Based on the above scenarios, the flow parameters are assigned to every individual node. At the same time, the thermal power denotes nozzles that are jetting.
denotes the adjacent nozzles that are used as thermal actuators. denotes the nozzles that are not used for heating. (a) Controlled scenario: left half of NP left is used for jetting. The maximum absolute difference of temperature between two adjacent nozzles is 0.1988 • C. (b) Uncontrolled scenario: left half of NP left is used for jetting. The maximum absolute difference of temperature between two adjacent nozzles is 1.2292 • C. dissipated by the individual jetting nozzle is used as known disturbances.
The model of every individual node is considered in discrete time by choosing T = {kt d | k ∈ N} with t d = 0.01 s. Euler's approach is used for time discretization due to its sparse and structure preserving implementation [28, p. 4].
The thermofluidic model is defined as a graph following the definitions (1)-(3) in Section III. The equivalent representations P II and P III are determined by eliminating the interconnection and print-job signals.
The in situ sensing-actuation scheme, as discussed in Section VI-A, is implemented by allocating the nozzles based on their four modes operations. The set S p (t) determines the set of all jetting nozzles. The rest of nozzles, except for the three nozzles that are located at the furthest distance from the jetting nozzles, are used for heating. This determines S h (t). All the nozzles that do not belong to S p (t) are used as soft sensors to estimate the change in liquid temperature. With the allocated nozzles, the matrices S h (t), S p (t), S s (t), and S u (t) are built. The model for control (16) is built subsequently. The unmeasured states are automatically replaced by the corresponding state updates from the model. (16), the prediction model is built over the time horizon t ∈ T k N , where T k N := {kt d | k ∈ N [k,N +k] }. During designing MPC, the tradeoff between computational expense and accuracy is carefully studied to select the horizon length. By performing various simulation studies, the horizon length is finally chosen as N = 8. In (20), the weights Q and R i|k are chosen as diagonal matrices. Here, the diagonal entries in R i|k are chosen significantly higher than that of Q to strictly penalize the deviation of inputs from its reference values. Moreover, entries of R i|k are chosen such that the heating nozzles adjacent to the jetting nozzles have more input power.

4) Solving MPC Scheme: Using
The optimization problem in (22) is solved by the freely available mpcqpsolver in MATLAB by using the interior point method. Once the optimization yields the optimal control inputs, they are applied to the heating piezoelectric actuators by means of nonjetting voltage pulses. The MPC then repeats the same procedure over the entire bitmap iteratively.

5) Results:
Once the MPC is applied, the liquid temperatures of all 48 nozzles are shown in Figs. [15][16][17]. It demonstrates that the performance specifications are met while satisfying the constraints. Moreover, only the adjacent heating nozzles are used to compensate for the temperature inconsistencies among jetting and nonjetting nozzles.
Due to the modularity of the digital twin, upscaling (or changing) the number of nozzles does not require rebuilding the entire model. The digital twin and the control software are automated to build an upscaled model with a user-defined number of nozzles, implement in situ sensing-actuation scheme using the bitmap, and visualize the results of the closed-loop system once the MPC is applied. To demonstrate that, the same printhead, shown in Fig. 1(a), is equipped with 160 nozzles that are equally divided over NP left and NP right . In each side, there are 80 nozzles that are equally distributed over four rows. In Fig. 18, the result of MPC applied configuration is shown when half the nozzles in NP left are used for jetting.

Remark 2 (On Computational Complexity):
With an increasing number of the nozzles, evidently, the computational complexity of solving MPC increases. At a specific iteration of k, the number of decision variables is related to the number of states, the number of control inputs, and the number of horizons. Let these numbers be n x , n k h , and N. In contrast to solving (20) and (21) that has a complexity of O(N 3 (n x +n k h ) 3 ), the dense LCQP (22) has a lower complexity O(N 3 n k h 3 ) [27]. This is achieved by eliminating the states from the decision variables with substitution of equality constraints. However, in this more condensed formulation, the sparsity of matrices is partially lost [29]. Yet, the dense LCQP (22) is to be preferred from a computational point of view as the dimension of decision variables involves only the number of nonjetting heating actuators. This is typically small in number.
For a few number of nozzles (below 30), it is possible to implement the controller during a printing operation. However, considering the large number of nozzles in an industrial printhead, direct implementation of the presented MPC may not be possible. A practical solution would compute the MPC controller gains offline (on the basis of the prior knowledge on the bitmap and flow parameters) and then apply the MPC as a feedforward action (e.g., using lookup table). Moreover, there are future possibilities on distributing the computation and hardware architecture to implement such a large-scale MPC.

VII. CONCLUSION
In this article, a modular and flexible digital twin is presented for modeling and control of thermofluidic processes in a DoD inkjet printhead. In particular, no additional sensors or actuators are incorporated by developing an in situ sensingactuation-based control strategy that minimizes the liquid temperature gradients among individual nozzle. To this end, an experimentally validated graph-theoretic modeling framework is developed that is modular up to an arbitrary number of nozzles. It is demonstrated that this model is flexible, scalable, and versatile and that a number of equivalent input-state-output representations can be derived in a straightforward and explicit manner from the model, depending on the intended application. A control strategy is implemented without using additional sensors and without using additional actuators. Specifically, to circumvent this limitation, the piezoelectric elements at every individual nozzle serves three roles: 1) it is a jetting actuator for depositing liquid; 2) it is a soft sensor for estimating liquid temperature; and 3) it is a control actuator to diminish gradient in liquid temperature among nozzles. The developed controller maintains the fluctuation of liquid temperature among nozzles well below a range of ±0.3 • C.
The tradeoff studies related to different control architectures for reducing signal overload and faster MPC computation are not addressed in this article. Moreover, designing and scheduling nonjetting voltage pulses that are required to apply control input without forming droplets of liquid are not discussed in this article.