# New Coordination Software for Parameter Identification Applied to Thermal Models of an Actuator Strut

Authors: Burkhard Hensel; Steffen Schroeder; Klaus Kabitzsch
Journal: Journal of Computational Engineering
 (2017)

---

## Abstract

The practical worth of models of technical processes depends on their accuracy, that is, the difference between model outputs and real measurements. For minimizing these differences, process identification methods are used. In this article, coordination software for process identification is presented which has the unique feature that it allows the integration of models that have been created with external tools, for example, Matlab or Python scripts. There is no need to transform the models into another type of software format to use the common identification coordinator. The concept of the software is described and two examples for the coupling with external simulation software are given. Additionally, this article contains a detailed case study of the parameter identification of two models using that identification coordination software. This highlights the benefit of the new coordination software regarding similar work flow for different model types. The modeled physical subject is the thermal behavior of an actuator strut.


---

## Body

## 1. Introduction

Models of technical systems are used for solving many kinds of problems reaching from prediction over simulation, condition monitoring, and failure detection to controller design. However, the quality of the results for all these goals depends on the model quality, that is, the accuracy. The smaller the difference between the model and the real process is, the more accurate a prediction, control loop, and so forth can be.All model types consist of a modelstructure and a number ofparameters. Thestructure defines the underlying equations, assumptions, and the resolution (granularity, order) of the model. One model structure can often be used for many technical processes with the samequalitative behavior. Theparameters (e.g., coefficients of polynomials, differential, or difference equations) can be adjusted for getting the desiredquantitative results. The task of finding the parameter values which minimize the difference between real measurements and simulation results of the model (for the same inputs and other conditions) is calledsystem identification,process identification, orparameter identification.There are a lot of tools which can be used for parameter identification; an overview is given in [1]. Probably the most frequently used software is Matlab [2]. There are many toolboxes for Matlab which support process identification for a lot of different kinds of models. However, for models which havenot been implemented in Matlab, usually other identification tools appropriate to the other modeling software must be used. This often results in suboptimal identification results as, for example, due to time reasons, only simple tuning methodologies are used, often just manual adjustments in a trial-and-error manner.Recently, higher level process identification software has been presented, which neither providesnew process identification methods norreplaces existing identification software. Instead, it sets a coordination view on top of them [1]. Besides process identification, this software has its focus on virtual sensor design and can thus be used not only on standard PCs in offline mode, but also on embedded PCs in online mode. In the state at the publication of [1] it was not possible to use external models in that software. The software has now been extended to support different external model types (e.g., Matlab and Python scripts) additionally to “internal” model types (strictly speaking realized as Java plug-ins) like static equations, difference equations, or discrete time and continuous-time transfer functions. This allows using the identification software for already existing models that have been implemented in external software without the need to convert these models to new software’s format. Additionally, it is even possible to combine external models (cosimulation). For example, the parameters of a Matlab script model can in that way easily be optimized simultaneously with parameters of a Python script model.The article is structured as follows. In Section2, the basic system identification procedure, which is the background of the software, is presented. An overview about the high level identification software is given in Section 3 with the focus on the integration of external models. In Section 4, an application example (case study) is presented in detail. In that example, two different model types are used to model the same physical object. The example is the thermal modeling of an actuator strut that can be found in similar form in many machine tools. Although the models have been implemented in different software, the new higher level software has been used for both models to coordinate the process identification. The section compares both the quality of the models and the effort for creating them. Finally, conclusions are drawn (Section 5).The abbreviations and symbols used in this article are presented in Notations.


## 2. Identification Procedure

In this section, the identification procedure is presented, which is supported by the software. It is based on [1, 3, 4].Lety_k be a vector of M time series where k is the index of the sampling instant. Each element of the vector represents a time series for one physical quantity. More concretely, y_modelk is the vector of the outputs of the model and y_measurementk the vector of measured time series for the same physical quantities and under the same conditions and inputs. The difference between the model outputs and the corresponding measurements is mostly evaluated using the sum of squared error (SSE)(1)SSE_=SSE1SSE2⋯SSENT,(2)SSEi=∑j=1Nymodel,ik-ymeasurement,ik2,where N is the number of recorded time steps (sampling instants), SSE_ is the vector of squared errors, and i is the row index in the vectors y_modelk, y_measurementk, and SSE_. Another often used measure is the root-mean-square error (RMSE)(3)RMSE_=RMSE1RMSE2⋯RMSENT,RMSEi=SSEiN.In general, there are two main strategies of parameter identification.(1)
The first one is theanalytical computation of the optimal parameters for a given model structure according to given measurements. For that purpose, there are mathematical strategies like least squares method, instrumental variable method, maximum likelihood method, and so forth [5, 6]. These methods deliver the optimal parameters for a given optimization criterion (usually SSE_ with M=1) after a finite number of computational steps. However, this is only possible for certain model structures, for example, linear models, polynomials for static models, or linear difference equations for dynamic models, mostly models that are linear in their parameters. Most text books as well as publications on process identification focus on analytical methods [5, 6].(2)
On the other hand, there is the possibility of “training” models usingiterative optimization methods. In this case, the model is simulated again and again with varied parameters and the difference (SSE_) between the model output and the measurements is evaluated. There are a lot of optimization strategies like Monte Carlo optimization, simulated annealing, Levenberg-Marquardt, gradient-based search, and evolutionary algorithms. Most such algorithms take less assumptions on the model structure than analytical methods and can therefore be used for a wide variety of model types or even unknown (black-box) models. Therefore, they are the usual way of process identification for complex model types like finite element models [7, 8]. But also for simple models, where analytical methods exist, such optimization methods are sometimes used; a well-known example is the Matlab System Identification Toolbox [9]. However, these methods do not guarantee to deliver the globally optimal results (due to local optima), especially not in finite time. Both types of identification procedures are supported by the new software to give a maximum of freedom to the user. It is also possible to combine both strategies. For example, linear model parts can first be identified analytically, and more complex or even black-box parts can be identified via evolutionary algorithms, afterwards.


## 3. Software Overview

A first description of the developed software has been given in [1]. The coordinator manages a modular model where the subsystems are connected via “signals,” as scientifically known from “systems theory.”Most functionality is provided to that core viaplug-ins. There are currently three types of plug-ins:(1)
System-processing plug-ins provide model type specific functions, especially model simulation, parameter identification, and inversion.(2)
Reader plug-ins are responsible for getting measurement data into the software—for example, from files, data bases, or (in online mode) directly from sensors—that is used to compare the model outputs with the real modeled system.(3)
Visualization plug-ins are used to visualize simulation results, for example, graphically or as tables, or for writing output data into files. Besides the management of the modular model and system identification work flow the core contains standard optimization methods that are independent of the model type and can thus be used for all model types together.One of the reasons for the plug-in concept is the ability to extend the identification software by functionality of other software products using “wrapper plug-ins.” That concept has now been used to integrate models that have been implemented as Matlab or Python scripts. For the end user of the software, the following work flow results.(1)
Usually, models that have been written as scripts contain adjustable parameters asvariables of the software. If that is not the case, variables for uncertain parameters that are to be identified have to be created in the script to use the identification coordinator.(2)
After that, in the GUI (graphical user interface) of the identification coordinator, a system of type “Matlab script” or “Python script” is created. In the corresponding dialog the used script file name and the names of the parameters that have to be optimized are edited (Figure1). It is also possible to set parameter values manually from the software, which is helpful if the same script is used in different contexts.(3)
The needed other parts for process identification are specified. This contains the selection of an appropriate “reader” (e.g., file format, data base) for the measurement data that is to be used, as well as—if necessary—additional preprocessing steps or further model parts that are not part of the script.(4)
Since scripts can usually contain all thinkable types of models or combinations of them, they are treated as black-box in the current implementation. The Matlab and Python plug-ins provide therefore no model-specific optimization methods. Because of this, the window for iterative standard parameter optimization can be opened. The parameters that shall be optimized automatically are selected together with their expected bounds and—if needed for the optimization method that is chosen later—start values (Figure2).(5)
In the next step, the optimization criterion is chosen. This is a set of signal pairs (usually a pair of one measurement signalymeasurement and one simulation output signal ymodel) together with a weighting factor for each pair. Let w_ be the vector of weighting factors with size M; then the error(4)e=w_·SSE_is the evaluated criterion, withSSE_ taken from (1).(6)
Now, the identification method is chosen. Only model independent identification methods are provided by the core, for example, Monte Carlo simulation, simulated annealing, and an evolutionary algorithm.(7)
Finally, the parameter optimization can be started. The possible stop criteria are currently a defined number of iteration steps and a time limit. A quality limit (i.e., desired value ofe) could be added easily. The parameter optimization finishes when (at least)one of the specified criteria is reached.The coupling with Matlab is realized with the “Matlab Engine API for Java” [10] and the coupling with Python by integrating the open-source Jython interpreter [11]; see Figure 3. These two solutions represent also two of the three main possible ways of integrating external software into a Java plug-in: (1)
Use an API (application programmer’s interface) of the external software to directly communicate with the external software for setting parameter values, starting simulation, and exchanging signals.(2)
Provide aninterpreter or simulator for the external models as a plug-in.(3)
Start a simulation as abatch job, that is, via command line. The last type of coupling should only be chosen if the other two solutions are not possible or if the required libraries are too large for the application (e.g., on embedded systems), because starting a simulation as a batch job is usually slower, especially if the start of the external simulation software needs significant time.Figure 1
Screenshot of the parameter selection for a script (variable names and values). The physical meaning of the parameters is explained in Section4.2.Figure 2
Screenshot of the form for selection of the parameters to be identified with their lower and upper limits. The physical meaning of the parameters and their limits is explained in Section4.2.Figure 3
Integration of external scripts in the identification coordinator.


## 4. Case Study: Actuator Strut

In this section, the concept of the software is explained using a practical application example, that is, the comparison of two model types for the same modeled object.The system to be modeled is an actuator strut. This example is taken from the collaborative research center Transregio 96, in which the software for parameter identification has been created. Actuator struts are typically used in machine tools, where for each axis (x, y, and z) one actuator strut is used to move either the work piece or the main spindle with the tool. Figure 4 shows the modeled actuator strut which has been separated from the overall machine for instrumentation and experimentation. It has only one relevant way of movement: the spindle can move out of and into the tube. Since one spindle end is fixed at the drive end (“bearing”), the spindle itself does not move, but the tube does (around the spindle). However, this has no influence on the thermal model as it is only a question of a coordinate transformation.Figure 4
Actuator strut with measurement points. “TS” = “temperature sensor.”From the research project’s perspective, the interesting quantity is the thermal “deformation” which is here only relevant in one dimension (axial), that is, the change of the length of the actuator strut due to temperature changes. This is interesting because thermal deformation in machine tools leads to less accuracy of production processes like milling or drilling. The thermal deformation is a function of the temperature field in the actuator strut. Therefore, the possibility of computing the temperature field is practically enough to estimate the change of the length that can be measured only with much higher costs than the temperature.Since the object is long and thin and only at one end fixed inx direction, only one direction (axial) is important for practical use. Assuming that for each point x on the axial direction the strut has the same temperature ϑx,t for each y and z, the change Δl of the length can be estimated from(5)Δl=∫0lαxΔϑxdx,where αx is the linear heat expansion coefficient for the substance at point x, Δϑx the change of the temperature at point x compared to its initial value, and l the initial length [12].The heat flow in a solid is described by the heat equation(6)c·∂2ϑ∂x2+∂2ϑ∂y2+∂2ϑ∂z2=∂ϑ∂t,where c is a constant (thermal diffusivity). However, this three-dimensional partial differential equation is too complex for practical simulation. At the same time, it is even not sufficient for describing the thermal behavior of the actuator strut, because of missing consideration of convection or heat radiation at surfaces. Therefore, several simplified model structures are available to approximate the behavior. Two of them will be used in the models described below. Theinput quantity of the models is the feed-rate (speed) and theoutputs are the temperatures at the measured points. The deformation (change of length) is left out of consideration in the models but could be added using (5).Starting from a thermal steady state, temperature curves (step responses) have been measured at several positions at the strut. For gaining these measurements, the actuator strut has been moved periodically with a constant feed-rate (velocity) for a duration between 80 and 170 minutes. Measurements have been taken for the feed-rates 20 m/min, 25 m/min, and 30 m/min. The sensor positions are shown in Figure4 and the measurement results in Figure 5.Figure 5
Measurements used for parameter identification.
(a)
v
=
20 m/min
(b)
v
=
25 m/min
(c)
v
=
30 m/minIn the following, it is shown how two models are set up using the data for a feed-rate of 20 m/min. The parameters are estimated using the new software. The quality of each model is evaluated by simulating the model for the feed-rates 20 m/min, 25 m/min, and 30 m/min and comparing the outputs with the appropriate measurements. The samples are nearly equidistant, but for technical reasons the sampling period depends on the feed-rate and is thus not equal for the three cases. For the casev= 20 m/min, the sampling period is 2.35 min (140.9 s), for v=25 m/min 2.48 min (148.7 s), and for v=30 m/min 2.64 min (158.3 s).


### 4.1. Model 1: Set of Transfer Functions

The measured temperature curves of the experiments shown in Figure5 are typical step responses of nonintegrating systems; that is, the influence of the speed on each temperature measurement point can be modeled either by a higher order dynamic system, a series of first-order systems, or (roughly) a first-order system with a time delay (FOPTD). The latter has the advantage that the parameters can be interpreted intuitively (but not physically) by comparing them with the step response. The transfer function of a FOPTD model is(7)Gis=ΘisVs=Ki1+Tise-sLi,where Ki is the proportional action coefficient, Ti the time constant (first-order lag), Li the time delay or dead time, and s the Laplace variable (complex frequency of the Laplace transform [13]). Θis is the Laplace transform of a (scalar) output (one temperature measurement point ϑis). Vs is the Laplace transform of the input (feed-rate vt). i is the index of the measurement position. The transfer function (7) is equivalent to the differential equation(8)ϑit+Tidϑitdt=Ki·vt-Li, where t is the (continuous) time.Even this simple model type has a relatively complex modular structure (see Figure6), because for each measurement point a separate FOPTD model has to be set up and the initial temperature of the measurements has to be subtracted. The latter is the case, because for the input (speed) equal to zero a linear model like an FOPTD model always delivers zero as output.Figure 6
Structure of the model consisting of FOPTD models.Since there are seven measurement points and each FOPTD model contains three parameters, the overall model has totally 21 parameters.The values of the parametersKi and Ti can be identified using the least squares method [5, 6, 14]. For that purpose, first the coefficients am1 and bm0 of a discrete difference equation(9)ϑik+am1·ϑik-1=bm0·vk-τ~ or the equivalent transfer function(10)Gdz-1=z-τ~bm01+am1z-1=z-τ~Ki·1-m11-m1z-1with the discrete time delay τ~=Li/Ts and(11)m1=exp⁡-TsTi,am1=-m1,bm0=Ki·1-m1 are computed. z is the variable of the z-Transform [13] and Ts the sampling period. Then, the continuous-time parameters Ki and Ti can be computed from them using(12)K^i=Gd1=bm01+am1,T^i=-Tsln⁡m1=-Tsln⁡-am1.The discrete time delayτ~ (and thus also the continuous-time time delay Li) cannot be computed using the least squares method. Therefore, for each integer value of τ~ in the interval 0,20 the model parameters have been computed (automatically) and the one with the smallest mean quadratic error has been used.This model type is not provided via external software as it is relatively simple. The model structure, the appropriate version of the least squares algorithm, and the loop of time delays have been directly programmed as a Java plug-in without a link to an external software product.Table1 shows the identified values of the parameters for each measurement point for the feed-rate 20 m/min.Table 1
FOPTD models for temperature sensors. Optimized forv=20 m/min (sampling period 2.35 min), other cases as validation without new parameter optimization.
Number
Sensor
DelayL in min
Prop. act. coeff.K in Kmin/m
Time constantT in min
RMSE 20 m/minin K
RMSE 25 m/minin K
RMSE 30 m/minin K
8
Spindle
0
1.33
22.39
0.66
2.49
7.96
5
Ball screw nut
0
1.24
28.25
0.74
1.78
7.60
3
Drive
0
1.20
39.06
0.45
5.59
7.05
4
Bearing
0
1.05
42.52
0.40
2.48
5.20
2
Belt housing
4.70
0.85
57.53
0.51
1.38
3.83
6
Tube middle
2.35
0.52
41.99
0.48
0.42
2.50
7
Tube end
11.74
0.15
93.04
0.08
0.19
0.56
Mean of all
—
—
—
0.47
2.05
4.96The validation cases show much larger values ofRMSE than the optimized case. This is a hint that the true modeled object is nonlinear. In the case v=25 m/min, the biggest differences occur for the drive, because it had not fully cooled down before the experiment started. For v=30 m/min all values of RMSE become much larger. This can be explained with the nonlinearity, as according to the measured data the final temperature does not change significantly although the final temperature of the model output is proportional to v due to the superposition principle of linear models.Positions that are far away from any heat source show the largest time constants, the highest delay, and the smallest proportional action coefficient. The curve of the belt housing is typical for a higher order process due to the long way from the heat sources to the measurement point. Since the model is not of higher order, the time delay is in that case larger. Also the tube is reached relatively slowly by the heat, resulting in long time constants and a time delay.Since the delay can only be estimated as an integer multiple of the sampling period, further deviations occur. However, since the changes of the sampling period are below 13%, this influence can be neglected compared to large influence of the nonlinearity. Also the importance of the influence of the air temperature is minor.


### 4.2. Model 2: Electric Equivalent Circuit Diagram

#### 4.2.1. Model Overview

In this section, the actuator strut is modeled using an electric equivalent circuit diagram (based on electrothermomechanical analogies), shown in a simplified form in Figure7. In an electric equivalent circuit diagram, each relevant physical component (with index i) of the actuator strut is modeled by a node with a heat capacity Ci and a power loss Q˙i (heat source). Between each pair of nodes i,j there is thermal conductance gij. If two nodes have no direct contact to each other, the conductance is zero.Figure 7
Electric equivalent circuit diagram of the actuator strut. “Resistors” represent thermal conductance. The circles are simplified visualizations for more complex electric circuits consisting of a parallel circuit of a current source (heat loss) and a (heat) capacity between the node and “electric ground.” The size of a circle correlates with the size of the capacity, the size of an arrow correlates with the size of the current source, and the size of a resistor correlates with the size of conductance.Each nodei of the P nodes is modeled with the first-order ordinary differential equation(13)Ci·dϑitdt=Q˙i+∑j=1,j≠iPgij·ϑjt-ϑit.However, in practice, the heat flowsQ˙it, which are the inputs of the model (13), are not known. Instead, the feed-rate vt is known and it is assumed that the heat flows Q˙it are a function of vt. In the simplest case, a proportional dependency can be assumed; that is, (14)Q˙it=qi·vtwith a proportionality factor qi (in J/m).


#### 4.2.2. Parameters of the Model

The model contains three types of parameters: a heat capacityCi and a velocity-dependent heat loss qi for each node and heat conductance gij between each pair of connected nodes. Nodes that are not connected have a heat conductance gij of 0. The heat capacity of the air can be treated as infinity.Sincegii is not needed and gij=gji, it is enough to estimate the heat conductance gij with i<j, that is, P-1·P/2 heat conductance parameters. Together with the P heat capacities and P heat losses, the total number of parameters can be computed to 2·P+P-1·P/2. The model of Figure 7 contains 19 nodes, so the number of parameters is 209.


#### 4.2.3. Parameter Reduction

The whole model can be represented by a system of first-order differential equations according to (13). Therefore, it is possible to create an equivalent state-space model and use analytical methods for parameter identification [15, 16]. However, the analytical process identification methods for state-space models are quite complex, because not only are there many parameters to be estimated (here 209), but also the state has to be estimated, here the temperatures of all model nodes—also these that are not measured at the real machine. Additionally, in the present example the model has been implemented as a Python script (see Figure 8) and the goal is to reuse this script and not to convert it into a state-space model of another type of software. In that way, this script serves as an example for the integration of arbitrary external scripts onto the identification coordinator. Because of that, only standard optimization methods are used. However, to identify 209 parameters simultaneously, especially with standard optimization methods, is practically hopeless because of the huge search space. Fortunately, the number of parameters to be identified can be heavily reduced using a priori knowledge. In particular, the 5 most uncertain parameters have been determined and only these are optimized by the optimization algorithm. All other parameters are set to values that have been estimated using physical knowledge. The strategy of finding the 5 most uncertain parameters and how the others have been estimated is described in Appendix.Figure 8
Systems and signals for electric equivalent circuit as a Python script.The remaining 5 parameters are “ratio parameters”; that is, they have no physical meaning but are factors that are multiplied with theoretically estimated physical parameters. The ratio parameterrm is related to the motor heat loss, rb to the heat loss of the bearing, rn to the heat loss of the ball screw nut, rsa to the heat conductance between the spindle and the air, and rta to the heat conductance between the tube and the air. See Appendix for details.The limits of the ratio parameters can be computed from the uncertainties shown in [4]; that is, the uncertainty is +200/−70% for the power loss of ball screw nut and bearing, +30/−20% for heat conductance with free convection and radiation, and ±25% for the motor power loss. Therefore(15)rm∈0.75,1.25,rb,rn∈0.6,2.0,rsa,rta∈0.8,1.3. In conclusion, the Python script has the five parameters rn, rb, rm, rsa, and rsb which are multiplied with the theoretically estimated model parameters according to Table 5.


#### 4.2.4. Evaluation

Each parameterri has been optimized in the given interval using Monte Carlo optimization. This optimization strategy is relatively slow but avoids sticking in local minima. The results are presented in Table 2. As for the FOLPD models, the parameters have been optimized for the data of v=20 m/min and the results of the unchanged model for v=25 m/min and v=30 m/min are given for comparison. The accuracy of the model is quite similar to the set of transfer functions. Especially regarding the validation cases (v=25 m/min and v=30 m/min) the differences are negligible. The reason is that both models are linear and thus suffer in an equal manner from the inability to model the nonlinear behavior of the physical process.Table 2
Optimal parameters of the equivalent circuit model for different optimization goals (optimized forv=20 m/min data). Optimization time 60 min per case and Monte Carlo optimization (about 7600–9200 iterations per case).
Optimized for sensor
r
n
r
b
r
s
a
r
t
a
r
m
SSE 20 m/minin K2
RMSE 20 m/minin K
RMSE 25 m/minin K
RMSE 30 m/minin K
Spindle
1.13
1.68
1.25
1.26
0.96
1.21
0.13
2.07
7.30
Ball screw nut
0.88
1.21
1.16
1.25
1.23
1.90
0.16
1.83
7.30
Drive
0.86
1.02
1.17
1.06
1.24
48.52
0.81
4.59
6.35
Bearing
1.74
0.70
1.13
1.12
0.92
6.10
0.29
2.40
4.93
Belt housing
1.58
0.76
1.00
1.03
0.77
6.79
0.30
2.03
4.13
Tube middle
0.96
0.99
1.07
0.94
1.06
0.74
0.10
0.29
2.43
Tube end
0.61
1.75
0.85
0.99
1.15
5.48
0.27
0.45
0.79
Mean of all
0.89
0.71
0.88
1.14
1.16
268.12
1.73
2.34
4.79It is evident that although for each simulation the same model structure of the overall machine has been used, the optimal parameters are very different depending on the quantity to be optimized. However, it is often possible that different settings produce the same output—which is one reason for this diversity. Others are given later in this section.As explained before, only 5 parameters have been optimized with optimization algorithms while in Section4.1 for the set of transfer functions 21 parameters have been optimized. Nevertheless, in the transfer function model only 3 parameters have an influence on one output signal while in the equivalent circuit model all 5 parameters have an effect on all output signals. These can be interpreted as degrees of freedom. Due to the five parameters, there is more freedom in the equivalent circuit model for adjusting the signal behavior at a given point as with the three parameters of the FOPTD model. Therefore, the accuracy of the equivalent circuit model is mostly higher than the accuracy of the FOPTD models when optimized for one signal to be as close as possible to the appropriate measured signal. However, there are also several shortcomings of this model type:(i)
The effort for making the model structure, implementing the script, computing the start values of the 209 parameters, and identifying the 5 most unknown parameters is much higher.(ii)
The time for optimization is much higher (60 minutes for the Monte Carlo algorithm with the numerous Python script calls versus less than 1 second for the least squares method with the FOPTD models) since iterative optimization methods have to be used instead of the least squares method.(iii)
If not only one signal but the average of theRMSEs of all 7 measurement points is optimized, then the model is much worse than the set of independent FOLPD models. One reason is that in the set of transfer functions 21 parameters have been optimized but for the equivalent circuit model only 5. Additionally, the transfer functions are not coupled and can thus be optimized independently from each other while in the equivalent circuit model each parameter has an influence on every output signal.(iv)
While one advantage of the equivalent circuit model is the physical interpretability of the parameters (capacities, heat conductivities, and heat losses), it is a wrong conclusion that the optimizer would bring the parameter values to the real physical parameters. Instead, it choosesany values that minimize the difference between measurements and simulation outputs, even if they are not physically correct but a result of the simplified model structure (a finite number of mass points) or unmodeled disturbances. The extreme differences between the optimal values of the parameters for the different optimization criteria (measurement points) confirm that. A further confirmation of that problem is the poor model quality when not only one signal but the average of all RMSEs is optimized. The conclusion of this problem is that the advantage of the equivalent circuit model—where the parameters are physically interpretable—is only working if the model structure is sufficiently correct and all parameters can be computed without parameter identification techniques. As soon as parameter identification is used, the identified parameters have no better physical meaning than the parameters of the phenomenological FOPTD model.Nevertheless, the choice of a model type always depends on the concrete application and the questions to be answered.


## 4.1. Model 1: Set of Transfer Functions

The measured temperature curves of the experiments shown in Figure5 are typical step responses of nonintegrating systems; that is, the influence of the speed on each temperature measurement point can be modeled either by a higher order dynamic system, a series of first-order systems, or (roughly) a first-order system with a time delay (FOPTD). The latter has the advantage that the parameters can be interpreted intuitively (but not physically) by comparing them with the step response. The transfer function of a FOPTD model is(7)Gis=ΘisVs=Ki1+Tise-sLi,where Ki is the proportional action coefficient, Ti the time constant (first-order lag), Li the time delay or dead time, and s the Laplace variable (complex frequency of the Laplace transform [13]). Θis is the Laplace transform of a (scalar) output (one temperature measurement point ϑis). Vs is the Laplace transform of the input (feed-rate vt). i is the index of the measurement position. The transfer function (7) is equivalent to the differential equation(8)ϑit+Tidϑitdt=Ki·vt-Li, where t is the (continuous) time.Even this simple model type has a relatively complex modular structure (see Figure6), because for each measurement point a separate FOPTD model has to be set up and the initial temperature of the measurements has to be subtracted. The latter is the case, because for the input (speed) equal to zero a linear model like an FOPTD model always delivers zero as output.Figure 6
Structure of the model consisting of FOPTD models.Since there are seven measurement points and each FOPTD model contains three parameters, the overall model has totally 21 parameters.The values of the parametersKi and Ti can be identified using the least squares method [5, 6, 14]. For that purpose, first the coefficients am1 and bm0 of a discrete difference equation(9)ϑik+am1·ϑik-1=bm0·vk-τ~ or the equivalent transfer function(10)Gdz-1=z-τ~bm01+am1z-1=z-τ~Ki·1-m11-m1z-1with the discrete time delay τ~=Li/Ts and(11)m1=exp⁡-TsTi,am1=-m1,bm0=Ki·1-m1 are computed. z is the variable of the z-Transform [13] and Ts the sampling period. Then, the continuous-time parameters Ki and Ti can be computed from them using(12)K^i=Gd1=bm01+am1,T^i=-Tsln⁡m1=-Tsln⁡-am1.The discrete time delayτ~ (and thus also the continuous-time time delay Li) cannot be computed using the least squares method. Therefore, for each integer value of τ~ in the interval 0,20 the model parameters have been computed (automatically) and the one with the smallest mean quadratic error has been used.This model type is not provided via external software as it is relatively simple. The model structure, the appropriate version of the least squares algorithm, and the loop of time delays have been directly programmed as a Java plug-in without a link to an external software product.Table1 shows the identified values of the parameters for each measurement point for the feed-rate 20 m/min.Table 1
FOPTD models for temperature sensors. Optimized forv=20 m/min (sampling period 2.35 min), other cases as validation without new parameter optimization.
Number
Sensor
DelayL in min
Prop. act. coeff.K in Kmin/m
Time constantT in min
RMSE 20 m/minin K
RMSE 25 m/minin K
RMSE 30 m/minin K
8
Spindle
0
1.33
22.39
0.66
2.49
7.96
5
Ball screw nut
0
1.24
28.25
0.74
1.78
7.60
3
Drive
0
1.20
39.06
0.45
5.59
7.05
4
Bearing
0
1.05
42.52
0.40
2.48
5.20
2
Belt housing
4.70
0.85
57.53
0.51
1.38
3.83
6
Tube middle
2.35
0.52
41.99
0.48
0.42
2.50
7
Tube end
11.74
0.15
93.04
0.08
0.19
0.56
Mean of all
—
—
—
0.47
2.05
4.96The validation cases show much larger values ofRMSE than the optimized case. This is a hint that the true modeled object is nonlinear. In the case v=25 m/min, the biggest differences occur for the drive, because it had not fully cooled down before the experiment started. For v=30 m/min all values of RMSE become much larger. This can be explained with the nonlinearity, as according to the measured data the final temperature does not change significantly although the final temperature of the model output is proportional to v due to the superposition principle of linear models.Positions that are far away from any heat source show the largest time constants, the highest delay, and the smallest proportional action coefficient. The curve of the belt housing is typical for a higher order process due to the long way from the heat sources to the measurement point. Since the model is not of higher order, the time delay is in that case larger. Also the tube is reached relatively slowly by the heat, resulting in long time constants and a time delay.Since the delay can only be estimated as an integer multiple of the sampling period, further deviations occur. However, since the changes of the sampling period are below 13%, this influence can be neglected compared to large influence of the nonlinearity. Also the importance of the influence of the air temperature is minor.


## 4.2. Model 2: Electric Equivalent Circuit Diagram

### 4.2.1. Model Overview

In this section, the actuator strut is modeled using an electric equivalent circuit diagram (based on electrothermomechanical analogies), shown in a simplified form in Figure7. In an electric equivalent circuit diagram, each relevant physical component (with index i) of the actuator strut is modeled by a node with a heat capacity Ci and a power loss Q˙i (heat source). Between each pair of nodes i,j there is thermal conductance gij. If two nodes have no direct contact to each other, the conductance is zero.Figure 7
Electric equivalent circuit diagram of the actuator strut. “Resistors” represent thermal conductance. The circles are simplified visualizations for more complex electric circuits consisting of a parallel circuit of a current source (heat loss) and a (heat) capacity between the node and “electric ground.” The size of a circle correlates with the size of the capacity, the size of an arrow correlates with the size of the current source, and the size of a resistor correlates with the size of conductance.Each nodei of the P nodes is modeled with the first-order ordinary differential equation(13)Ci·dϑitdt=Q˙i+∑j=1,j≠iPgij·ϑjt-ϑit.However, in practice, the heat flowsQ˙it, which are the inputs of the model (13), are not known. Instead, the feed-rate vt is known and it is assumed that the heat flows Q˙it are a function of vt. In the simplest case, a proportional dependency can be assumed; that is, (14)Q˙it=qi·vtwith a proportionality factor qi (in J/m).


### 4.2.2. Parameters of the Model

The model contains three types of parameters: a heat capacityCi and a velocity-dependent heat loss qi for each node and heat conductance gij between each pair of connected nodes. Nodes that are not connected have a heat conductance gij of 0. The heat capacity of the air can be treated as infinity.Sincegii is not needed and gij=gji, it is enough to estimate the heat conductance gij with i<j, that is, P-1·P/2 heat conductance parameters. Together with the P heat capacities and P heat losses, the total number of parameters can be computed to 2·P+P-1·P/2. The model of Figure 7 contains 19 nodes, so the number of parameters is 209.


### 4.2.3. Parameter Reduction

The whole model can be represented by a system of first-order differential equations according to (13). Therefore, it is possible to create an equivalent state-space model and use analytical methods for parameter identification [15, 16]. However, the analytical process identification methods for state-space models are quite complex, because not only are there many parameters to be estimated (here 209), but also the state has to be estimated, here the temperatures of all model nodes—also these that are not measured at the real machine. Additionally, in the present example the model has been implemented as a Python script (see Figure 8) and the goal is to reuse this script and not to convert it into a state-space model of another type of software. In that way, this script serves as an example for the integration of arbitrary external scripts onto the identification coordinator. Because of that, only standard optimization methods are used. However, to identify 209 parameters simultaneously, especially with standard optimization methods, is practically hopeless because of the huge search space. Fortunately, the number of parameters to be identified can be heavily reduced using a priori knowledge. In particular, the 5 most uncertain parameters have been determined and only these are optimized by the optimization algorithm. All other parameters are set to values that have been estimated using physical knowledge. The strategy of finding the 5 most uncertain parameters and how the others have been estimated is described in Appendix.Figure 8
Systems and signals for electric equivalent circuit as a Python script.The remaining 5 parameters are “ratio parameters”; that is, they have no physical meaning but are factors that are multiplied with theoretically estimated physical parameters. The ratio parameterrm is related to the motor heat loss, rb to the heat loss of the bearing, rn to the heat loss of the ball screw nut, rsa to the heat conductance between the spindle and the air, and rta to the heat conductance between the tube and the air. See Appendix for details.The limits of the ratio parameters can be computed from the uncertainties shown in [4]; that is, the uncertainty is +200/−70% for the power loss of ball screw nut and bearing, +30/−20% for heat conductance with free convection and radiation, and ±25% for the motor power loss. Therefore(15)rm∈0.75,1.25,rb,rn∈0.6,2.0,rsa,rta∈0.8,1.3. In conclusion, the Python script has the five parameters rn, rb, rm, rsa, and rsb which are multiplied with the theoretically estimated model parameters according to Table 5.


### 4.2.4. Evaluation

Each parameterri has been optimized in the given interval using Monte Carlo optimization. This optimization strategy is relatively slow but avoids sticking in local minima. The results are presented in Table 2. As for the FOLPD models, the parameters have been optimized for the data of v=20 m/min and the results of the unchanged model for v=25 m/min and v=30 m/min are given for comparison. The accuracy of the model is quite similar to the set of transfer functions. Especially regarding the validation cases (v=25 m/min and v=30 m/min) the differences are negligible. The reason is that both models are linear and thus suffer in an equal manner from the inability to model the nonlinear behavior of the physical process.Table 2
Optimal parameters of the equivalent circuit model for different optimization goals (optimized forv=20 m/min data). Optimization time 60 min per case and Monte Carlo optimization (about 7600–9200 iterations per case).
Optimized for sensor
r
n
r
b
r
s
a
r
t
a
r
m
SSE 20 m/minin K2
RMSE 20 m/minin K
RMSE 25 m/minin K
RMSE 30 m/minin K
Spindle
1.13
1.68
1.25
1.26
0.96
1.21
0.13
2.07
7.30
Ball screw nut
0.88
1.21
1.16
1.25
1.23
1.90
0.16
1.83
7.30
Drive
0.86
1.02
1.17
1.06
1.24
48.52
0.81
4.59
6.35
Bearing
1.74
0.70
1.13
1.12
0.92
6.10
0.29
2.40
4.93
Belt housing
1.58
0.76
1.00
1.03
0.77
6.79
0.30
2.03
4.13
Tube middle
0.96
0.99
1.07
0.94
1.06
0.74
0.10
0.29
2.43
Tube end
0.61
1.75
0.85
0.99
1.15
5.48
0.27
0.45
0.79
Mean of all
0.89
0.71
0.88
1.14
1.16
268.12
1.73
2.34
4.79It is evident that although for each simulation the same model structure of the overall machine has been used, the optimal parameters are very different depending on the quantity to be optimized. However, it is often possible that different settings produce the same output—which is one reason for this diversity. Others are given later in this section.As explained before, only 5 parameters have been optimized with optimization algorithms while in Section4.1 for the set of transfer functions 21 parameters have been optimized. Nevertheless, in the transfer function model only 3 parameters have an influence on one output signal while in the equivalent circuit model all 5 parameters have an effect on all output signals. These can be interpreted as degrees of freedom. Due to the five parameters, there is more freedom in the equivalent circuit model for adjusting the signal behavior at a given point as with the three parameters of the FOPTD model. Therefore, the accuracy of the equivalent circuit model is mostly higher than the accuracy of the FOPTD models when optimized for one signal to be as close as possible to the appropriate measured signal. However, there are also several shortcomings of this model type:(i)
The effort for making the model structure, implementing the script, computing the start values of the 209 parameters, and identifying the 5 most unknown parameters is much higher.(ii)
The time for optimization is much higher (60 minutes for the Monte Carlo algorithm with the numerous Python script calls versus less than 1 second for the least squares method with the FOPTD models) since iterative optimization methods have to be used instead of the least squares method.(iii)
If not only one signal but the average of theRMSEs of all 7 measurement points is optimized, then the model is much worse than the set of independent FOLPD models. One reason is that in the set of transfer functions 21 parameters have been optimized but for the equivalent circuit model only 5. Additionally, the transfer functions are not coupled and can thus be optimized independently from each other while in the equivalent circuit model each parameter has an influence on every output signal.(iv)
While one advantage of the equivalent circuit model is the physical interpretability of the parameters (capacities, heat conductivities, and heat losses), it is a wrong conclusion that the optimizer would bring the parameter values to the real physical parameters. Instead, it choosesany values that minimize the difference between measurements and simulation outputs, even if they are not physically correct but a result of the simplified model structure (a finite number of mass points) or unmodeled disturbances. The extreme differences between the optimal values of the parameters for the different optimization criteria (measurement points) confirm that. A further confirmation of that problem is the poor model quality when not only one signal but the average of all RMSEs is optimized. The conclusion of this problem is that the advantage of the equivalent circuit model—where the parameters are physically interpretable—is only working if the model structure is sufficiently correct and all parameters can be computed without parameter identification techniques. As soon as parameter identification is used, the identified parameters have no better physical meaning than the parameters of the phenomenological FOPTD model.Nevertheless, the choice of a model type always depends on the concrete application and the questions to be answered.


## 4.2.1. Model Overview

In this section, the actuator strut is modeled using an electric equivalent circuit diagram (based on electrothermomechanical analogies), shown in a simplified form in Figure7. In an electric equivalent circuit diagram, each relevant physical component (with index i) of the actuator strut is modeled by a node with a heat capacity Ci and a power loss Q˙i (heat source). Between each pair of nodes i,j there is thermal conductance gij. If two nodes have no direct contact to each other, the conductance is zero.Figure 7
Electric equivalent circuit diagram of the actuator strut. “Resistors” represent thermal conductance. The circles are simplified visualizations for more complex electric circuits consisting of a parallel circuit of a current source (heat loss) and a (heat) capacity between the node and “electric ground.” The size of a circle correlates with the size of the capacity, the size of an arrow correlates with the size of the current source, and the size of a resistor correlates with the size of conductance.Each nodei of the P nodes is modeled with the first-order ordinary differential equation(13)Ci·dϑitdt=Q˙i+∑j=1,j≠iPgij·ϑjt-ϑit.However, in practice, the heat flowsQ˙it, which are the inputs of the model (13), are not known. Instead, the feed-rate vt is known and it is assumed that the heat flows Q˙it are a function of vt. In the simplest case, a proportional dependency can be assumed; that is, (14)Q˙it=qi·vtwith a proportionality factor qi (in J/m).


## 4.2.2. Parameters of the Model

The model contains three types of parameters: a heat capacityCi and a velocity-dependent heat loss qi for each node and heat conductance gij between each pair of connected nodes. Nodes that are not connected have a heat conductance gij of 0. The heat capacity of the air can be treated as infinity.Sincegii is not needed and gij=gji, it is enough to estimate the heat conductance gij with i<j, that is, P-1·P/2 heat conductance parameters. Together with the P heat capacities and P heat losses, the total number of parameters can be computed to 2·P+P-1·P/2. The model of Figure 7 contains 19 nodes, so the number of parameters is 209.


## 4.2.3. Parameter Reduction

The whole model can be represented by a system of first-order differential equations according to (13). Therefore, it is possible to create an equivalent state-space model and use analytical methods for parameter identification [15, 16]. However, the analytical process identification methods for state-space models are quite complex, because not only are there many parameters to be estimated (here 209), but also the state has to be estimated, here the temperatures of all model nodes—also these that are not measured at the real machine. Additionally, in the present example the model has been implemented as a Python script (see Figure 8) and the goal is to reuse this script and not to convert it into a state-space model of another type of software. In that way, this script serves as an example for the integration of arbitrary external scripts onto the identification coordinator. Because of that, only standard optimization methods are used. However, to identify 209 parameters simultaneously, especially with standard optimization methods, is practically hopeless because of the huge search space. Fortunately, the number of parameters to be identified can be heavily reduced using a priori knowledge. In particular, the 5 most uncertain parameters have been determined and only these are optimized by the optimization algorithm. All other parameters are set to values that have been estimated using physical knowledge. The strategy of finding the 5 most uncertain parameters and how the others have been estimated is described in Appendix.Figure 8
Systems and signals for electric equivalent circuit as a Python script.The remaining 5 parameters are “ratio parameters”; that is, they have no physical meaning but are factors that are multiplied with theoretically estimated physical parameters. The ratio parameterrm is related to the motor heat loss, rb to the heat loss of the bearing, rn to the heat loss of the ball screw nut, rsa to the heat conductance between the spindle and the air, and rta to the heat conductance between the tube and the air. See Appendix for details.The limits of the ratio parameters can be computed from the uncertainties shown in [4]; that is, the uncertainty is +200/−70% for the power loss of ball screw nut and bearing, +30/−20% for heat conductance with free convection and radiation, and ±25% for the motor power loss. Therefore(15)rm∈0.75,1.25,rb,rn∈0.6,2.0,rsa,rta∈0.8,1.3. In conclusion, the Python script has the five parameters rn, rb, rm, rsa, and rsb which are multiplied with the theoretically estimated model parameters according to Table 5.


## 4.2.4. Evaluation

Each parameterri has been optimized in the given interval using Monte Carlo optimization. This optimization strategy is relatively slow but avoids sticking in local minima. The results are presented in Table 2. As for the FOLPD models, the parameters have been optimized for the data of v=20 m/min and the results of the unchanged model for v=25 m/min and v=30 m/min are given for comparison. The accuracy of the model is quite similar to the set of transfer functions. Especially regarding the validation cases (v=25 m/min and v=30 m/min) the differences are negligible. The reason is that both models are linear and thus suffer in an equal manner from the inability to model the nonlinear behavior of the physical process.Table 2
Optimal parameters of the equivalent circuit model for different optimization goals (optimized forv=20 m/min data). Optimization time 60 min per case and Monte Carlo optimization (about 7600–9200 iterations per case).
Optimized for sensor
r
n
r
b
r
s
a
r
t
a
r
m
SSE 20 m/minin K2
RMSE 20 m/minin K
RMSE 25 m/minin K
RMSE 30 m/minin K
Spindle
1.13
1.68
1.25
1.26
0.96
1.21
0.13
2.07
7.30
Ball screw nut
0.88
1.21
1.16
1.25
1.23
1.90
0.16
1.83
7.30
Drive
0.86
1.02
1.17
1.06
1.24
48.52
0.81
4.59
6.35
Bearing
1.74
0.70
1.13
1.12
0.92
6.10
0.29
2.40
4.93
Belt housing
1.58
0.76
1.00
1.03
0.77
6.79
0.30
2.03
4.13
Tube middle
0.96
0.99
1.07
0.94
1.06
0.74
0.10
0.29
2.43
Tube end
0.61
1.75
0.85
0.99
1.15
5.48
0.27
0.45
0.79
Mean of all
0.89
0.71
0.88
1.14
1.16
268.12
1.73
2.34
4.79It is evident that although for each simulation the same model structure of the overall machine has been used, the optimal parameters are very different depending on the quantity to be optimized. However, it is often possible that different settings produce the same output—which is one reason for this diversity. Others are given later in this section.As explained before, only 5 parameters have been optimized with optimization algorithms while in Section4.1 for the set of transfer functions 21 parameters have been optimized. Nevertheless, in the transfer function model only 3 parameters have an influence on one output signal while in the equivalent circuit model all 5 parameters have an effect on all output signals. These can be interpreted as degrees of freedom. Due to the five parameters, there is more freedom in the equivalent circuit model for adjusting the signal behavior at a given point as with the three parameters of the FOPTD model. Therefore, the accuracy of the equivalent circuit model is mostly higher than the accuracy of the FOPTD models when optimized for one signal to be as close as possible to the appropriate measured signal. However, there are also several shortcomings of this model type:(i)
The effort for making the model structure, implementing the script, computing the start values of the 209 parameters, and identifying the 5 most unknown parameters is much higher.(ii)
The time for optimization is much higher (60 minutes for the Monte Carlo algorithm with the numerous Python script calls versus less than 1 second for the least squares method with the FOPTD models) since iterative optimization methods have to be used instead of the least squares method.(iii)
If not only one signal but the average of theRMSEs of all 7 measurement points is optimized, then the model is much worse than the set of independent FOLPD models. One reason is that in the set of transfer functions 21 parameters have been optimized but for the equivalent circuit model only 5. Additionally, the transfer functions are not coupled and can thus be optimized independently from each other while in the equivalent circuit model each parameter has an influence on every output signal.(iv)
While one advantage of the equivalent circuit model is the physical interpretability of the parameters (capacities, heat conductivities, and heat losses), it is a wrong conclusion that the optimizer would bring the parameter values to the real physical parameters. Instead, it choosesany values that minimize the difference between measurements and simulation outputs, even if they are not physically correct but a result of the simplified model structure (a finite number of mass points) or unmodeled disturbances. The extreme differences between the optimal values of the parameters for the different optimization criteria (measurement points) confirm that. A further confirmation of that problem is the poor model quality when not only one signal but the average of all RMSEs is optimized. The conclusion of this problem is that the advantage of the equivalent circuit model—where the parameters are physically interpretable—is only working if the model structure is sufficiently correct and all parameters can be computed without parameter identification techniques. As soon as parameter identification is used, the identified parameters have no better physical meaning than the parameters of the phenomenological FOPTD model.Nevertheless, the choice of a model type always depends on the concrete application and the questions to be answered.


## 5. Conclusions

There are two main goals of this article. The first is the presentation of a higher level software for parameter identification that allows now the inclusion of models that have been realized with external software without the need for remodeling in the new software. In particular, this has been realized for Matlab and Python scripts, demonstrating the generic concept. Due to the plug-in-based architecture, it is possible to extend the software for coupling a lot of other simulation software without changing the core of the identification coordinator. Additionally, this allows reusing established simulation and modeling software with all their specific features.The second goal of this article is the detailed presentation of a case study where two models of an actuator strut are compared with respect to parameter identification. A set of transfer functions has been compared to an equivalent circuit model. Both have their advantages and disadvantages regarding reachable accuracy, modeling effort, identification effort, and physical interpretability. Although with an equivalent circuit model a higher accuracy can be reached for single output signals, this results in less physical interpretability and less quality when all output signals are optimized. Both models have been identified using the new identification coordinator software such that the reading of measurement data and signal comparison “infrastructure” could be reused, saving time for parameter identification.



---

*Source: 3169785-2017-10-10.xml*