On Using Quasi-Newton Algorithms of the Broyden Class for Model-to-Test Correlation

ABSTRACT:

The correlation of a model with test results is a common task in engineering. Often genetic algorithms or adaptive particle swarm algorithms are used for this task. In this paper another approach is presented using two quasi-Newton algorithms of the class defined by Broyden. A study was conducted with thermal models showing the performance of this approach. Comparing the results to other studies reveals that the approach reduces the number of iterations by several orders of magnitude; typical calculation times for model correlation times are reduced from the order of weeks and months to the order of hours and days.

KEYWORDS:
Thermal vacuum test; Thermal analysis; Correlation; Broyden; ESATAN; Thermica

INTRODUCTION

In every calculation, there are always discrepancies between the predictions of a mathematical model and their physical measurements. Refinement of the models can reduce the discrepancies, however, the exact values of many parameters are often unknown. Without perfect knowledge of the values of these parameters, a model refinement will never lead exactly to the measured results.

Determining the exact values of all parameters individually is too costly or sometimes even impossible. Testing is often only possible at the system level, and extracting parameter values from the test results is a difficult task which is done by correlating the mathematical model with the measurement data. This means that the model parameters are changed via feedback from the measurement results, so that the discrepancies between the measurements and the models are minimized.

Many methods have been developed and analyzed to perform model-to-measurement correlation (Jouffroy, 2007; De Palo et al., 2011; Momayez et al., 2009; Harvatine and De Mauro, 1994; Roscher, 2006; van Zijl, 2013; WenLong et al., 2011; Mareschi et al., 2005). Most methods are based on stochastic optimization algorithms and often require several hundred iterations to converge. In this paper, an approach is presented using Broyden's class of methods (Broyden, 1965), which use considerably less iterations for the majority of cases.

THE THERMAL MODEL CORRELATION PROBLEM

The method will be presented with a typical thermal model correlation as used in the space industry, but the procedure can also be applied to many other models as used in other industries and areas.

Thermal models calculate temperatures at certain points depending on a set of parameters. Mathematically expressed, this means:

(1) $t_{\mathrm{mdl}}=F(p\mathrm{),}$

where:

F is the mathematical model with function F:Rk→Rm (k is the number of parameters, m is the number of results);

t_mdl is a vector containing the temperatures at the selected points; and

p is a vector containing the parameters of the model, for example thermal conductivity, thermal capacity, thermal emissivity, and solar absorptivity.

During the tests, the temperatures are measured at the same points as described in the model t_mdl. The aim of the correlation is to find the set of parameters p_corr for which the root sum square of the differences between the measured temperatures t_mes and the calculated temperatures F(p_corr) is minimal.

(2) $\mathrm{||}F(P_{\mathrm{corr}})t_{\mathrm{mes}}\mathrm{||}=\min \mathrm{\{||}F(p)-t_{\mathrm{mes}}\mathrm{|||}p\in P\}$

where P is the solution space for the parameter vector.

THE METHOD

In contrast to most correlation methods developed (Jouffroy, 2007; De Palo et al., 2011; Momayez et al., 2009; Harvatine and DeMauro, 1994; Roscher, 2006; van Zijl, 2013; WenLong et al., 2011; Mareschi et al., 2005), this method will not attempt to minimize the length of the vector ||F(p) - t_mes||, but searches a root of the equation system:

(3) $F(p_{\mathrm{corr}})-(t_{\mathrm{mes}})=0.$

The advantages of this method are that a whole vector of information, in this case temperature differences, can be evaluated instead of a single scalar and that a set of simple functions are used instead of a complex function. One of the best methods to solve this type of problem is the multidimensional Newton method, which applies the following formula iteratively:

(4) $p_{\mathit{n}+1}=p_{\mathit{n}}-\mathit{J}(p_{\mathit{n}}{)}^{-1}F(p_{\underset{─}{\mathit{n}}})$

where:

p_n is the vector of parameters;

J (p_n) is the Jacobian matrix at p_n

A limitation of this method is that the Jacobian matrix can often only be calculated numerically, with the consequence that, for each individual parameter in p_n, a separate calculation is needed. Fortunately, Broyden developed a class of methods to approximate the Jacobian matrix based on the results from the previous step. The methods are described in Broyden (1965) and therefore will not be explained here. Broyden's first method (the "good Broyden method") is chosen for the algorithm used in this paper as it performed better for Broyden's practical problems.

In addition, a self-developed method of the Broyden class is tested. This method updates each element b_n^ij of the approximated Jacobian Matrix B_n with the following formula:

(5) ${\mathit{b}}_{\mathit{n}+\mathit{l}}^{\mathit{ij}}=(1+{\mathit{k}}_{\mathit{n}}^{\mathit{i}}{\mathit{b}}_{\mathit{n}}^{\mathit{ij}}{\mathit{s}}_{\mathit{n}}^{\mathit{j}}){\mathit{b}}_{\mathit{n}}^{\mathit{ij}}$

where:

(6) ${\mathit{s}}_{\mathit{n}}:={\mathit{p}}_{\mathit{n}}-{\mathit{p}}_{\mathit{n}-1}$

(7) ${\mathit{y}}_{\mathit{ex},\mathit{n}}:={\mathit{B}}_{\mathit{n}}{\mathit{s}}_{\mathit{n}}$

(8) ${\mathit{y}}_{\mathit{is},\mathit{n}}:=\mathit{F}({\mathit{p}}_{\mathit{n}})-\mathit{F}({\mathit{p}}_{\mathit{n}-1})$

(9) ${\mathit{k}}^{\mathit{i}}=\frac{{\mathit{y}}_{\mathit{is},\mathit{n}}^{\mathit{i}}-{\mathit{y}}_{\mathit{ex},\mathit{n}}^{\mathit{i}}}{{\sum }_{\mathit{j}=1}^{\mathit{m}}({\mathit{b}}_{\mathit{n}}^{\mathit{ij}}{\mathit{s}}_{\mathit{n}}^{\mathit{j}}{)}^2}$

CONSTRAINTS FOR THE METHOD

Every optimization method has requirements which need to be fulfilled in order to be used effectively. Within this section, some of the more important requirements for practical model-to-test correlation, and the main practical consequences, are described. The requirements are important when choosing to use this method and also for the selection of the parameters and results to be used.

MONOTONY AND DIFFERENTIABILITY

An important restriction in the use of this method is that each result of the model is a differentiable monotone function of the parameters within the solution space.

Most thermal models are monotone and differentiable functions of the parameters, but there are a few exceptions such as controlled heating systems or when the parameter is, for example, the orbit position of a spacecraft. For small changes the model is neither monotone nor differentiable due to numerical noise.

If the function does not fulfill this requirement, within the solution space, the algorithm may be prevented from converging.

OBSERVABILITY OF THE PARAMETERS

Each of the parameters must have an effect on at least one result. This is not a specific issue to this particular method, but it is a general requirement for correlation. Although the requirement appears trivial, it is often the reason why the algorithm is unable to converge. For the presented method, a low degree of observability combined with numerical noise can lead to instability.

Often, the relevance of parameters can be deduced through an understanding of the model, but for complex systems it is sometimes difficult to estimate the relevance of an individual parameter. Fortunately, the Jacobian matrix addresses some of these limitations. If a column of the Jacobian matrix consists solely of zeros, the corresponding parameter cannot be observed. In practical applications there are often many values close to zero, mainly due to small effects or numerical noise. If the relevance of the parameter is smaller than the accuracy of the measurement data, the parameter can then be considered as not observable.

INFLUENCE OF THE PARAMETERS ON THE RESULTS

Each result has to be influenced by at least one parameter. This is also not a general requirement for correlation, but it prevents model correlation. For the presented method, numerical noise combined with a low level of influence on a result can lead to instability.

Again, the Jacobian matrix is a good tool for estimating the influence of a parameter on results. If a row of the Jacobian matrix has only zeros or very small values, the corresponding result cannot be influenced.

ACCURACY OF THE MEASUREMENT AND MODEL INACCURACIES

Due to measurement errors and model inaccuracies, it is unlikely that the result vector will converge towards zero. If the measurement error is larger than the effect of a parameter, then the relevant parameter may not converge to a reasonable result. Model inaccuracies may, in addition, make it impossible to correlate it within the parameter space.

CONSTRAINTS

In most instances, each parameter has certain limits defined by the physics of the model. For example, thermal conductivity will never be zero or below and never exceed a certain value. As the method defined by Broyden is for unconstrained problems, the following algorithm is used:

• Calculate p_n+1 without considering the constraints;

• If a parameter of p_n+1 is outside the defined boundaries, scale the change of the parameters so that all parameters are within the boundaries;

• If some parameters have already reached their constraints, calculate p_n+1 for the condition that these parameters are fixed at their constraints.

LOAD CASES ANALYZED

Two load cases were analyzed: a constructed simple model and a complex model of real piece of hardware. These models where build in Systema/THERMICA, a thermal analysis software specialized to simulate spacecraft (mostly compatible to ESATAN and similar to SINDA). The standard syntax for these models was used.

THE CONSTRUCTED SIMPLE MODEL

The thermal model is a steady state model consisting of four diffusive nodes and an environment node at 0°C. Each node has a surface of 0.1 m² with an emissivity of 1, which radiates to the environment. 10 W are dissipated in each node. Node 1 is connected to the environment node with a thermal conductance of 1 W/K. Nodes 1 to 4 are all connected to each other by 6 thermal linear links (GL) as shown in Fig. 1.

Figure 1
Thermal test model schematics.

In order for the temperatures at the nodes to correspond to the measured values, the parameter values of the thermal links need to be changed. Therefore, the equation system to be solved by the algorithm is:

(10) ${\mathit{T}}_{\mathrm{calculated}}(\mathit{GL})-{\mathit{T}}_{\mathrm{measured}}=0,$

where:

T_calculated is the function which represents the model;

GL is the vector containing GL₁ to GL₆; and

T_measured is the vector of the calculated temperatures for GL_i = (0.1 + i/100)W/K. This is assumed to be a measurement.

The initial estimates for all GL are 0.5 W/K.

Three load cases are derived from this problem:

• An underdetermined system: All links are variable. There are six parameters (links) for four values (temperatures);

• A determined system: Links 3 and 6 are constant (GL₃= 0.13W/K, GL₆ = 0.16W/K). There are four parameters (links) for four values (temperatures);

• An overdetermined system: Links 3, 5, and 6 are constant (GL₃= 0.13W/K, GL₅= 0.6W/K, GL₆= 0.16W/K). There are three parameters (links) for four values (temperatures), and GL₅ is chosen so that it is impossible to solve the system.

THE COMPLEX MODEL

The complex model is a model of real spacecraft hardware, which was tested in a thermo vacuum chamber. In total, 13 parameters where chosen; they consist of thermal boundary conditions, thermal conductances, and emissivities for the infrared spectrum.

Temperature readings were taken from 26 sensors at 6 test points over time, yielding 156 evaluated temperatures.

RESULTS

RESULTS OF THE UNDERDETERMINED SIMPLE MODEL

Figure 2 shows the development of the root sum square of all temperature differences for both analyzed algorithms. The values of some significant points shown in this diagram are listed in Table 1. The model quickly converged toward 0, with both algorithms reducing the RSS by approximately one order of magnitude every two to three iterations, until the numerical accuracy limit of 1E-5 was reached. After the 4^th iteration, the self-developed algorithm was faster, but after the 9^th iteration, the Broyden algorithm converged faster. In conclusion, the performance of both algorithms can be considered comparable for this case.

Figure 2
Root sum square of the temperature differences over iterations of the underdetermined simple model.

Table 1
Selected points of the algorithm convergence of the underdetermined simple model.

Algorithm	Iterations for Jacobian matrix generation	Nº of iterations after Jacobian matrix generation	RSS/K
Broyden	6	0 (initial)	4.261
		6	5.921E-03
		9	2.646E-05
Self-developed method	6	0 (initial)	4.261
		6	7.652E-04
		12	5.385E-05

RESULTS OF THE DETERMINED SIMPLE MODEL

In the same way as shown for the underdetermined simple model, Fig. 3 and Table 2 show the root sum square development for the determined simple model. The model converged slower than for the underdetermined problem. The Broyden algorithm needed an average of 11 steps per order of magnitude while the self-developed algorithm needed an average of only 5 of them.

Figure 3
Root sum square of the temperature differences over iterations of the determined simple model.

Table 2
Selected points of the algorithm convergence of the determined simple model.

Algorithm	Iterations for Jacobian matrix generation	Nº iterations after Jacobian matrix Nº iterations after Jacobian matrix generation	RSS/K
Broyden	4	0 (initial)	3.442
		6	1.024
		41	9.487E-05
Self-developed method	4	0 (initial)	3.442
		6	5.813E-02
		22	1.414E-05

RESULTS OF THE OVERDETERMINED SIMPLE MODEL

The overdetermined model converged, as expected, not to 0 but to a finite value (0.0375K), as is shown in Fig. 4 and Table 3. For this load case, the self-developed method converged considerably faster, with approximately 4 iterations per order of magnitude, while the Broyden algorithm used approximately 8 iterations per order of magnitude.

Figure 4
Root sum square of the temperature differences over iterations of the overdetermined simple model.

Table 3
Selected points of the algorithm convergence of the overdetermined simple model

Algorithm	Iterations for Jacobian matrix generation	Nº iterations after Jacobian matrix generation	RSS/K
Broyden	3	0 (initial)	3.442
		4	1.869
		15	3.760E-02
Self-developed method	3	0 (initial)	3.442
		4	1.869
		8	3.757E-02

RESULTS OF THE COMPLEX MODEL

Figure 5 and Table 4 show the convergence of the algorithm for the complex model in a different way to the way that was shown for the simple model. The complex model is an overdetermined system and it is not expected to converge to 0. Therefore, the difference between the RSS of each iteration and the best RSS reached is used. This difference is standardized to the initial value. The best RSS was reached by restarting the algorithm near the optimum.

Figure 5
relative root sum square of the temperature differences over iterations of the complex model.

Table 4
Selected points of the algorithm convergence of the complex model.

Algorithm	Iterations for Jacobian matrix generation	Nº iterations after Jacobian matrix generation	(RSSRSS min) / (RSSinitial- RSSmin)
Broyden	13	0 (initial)	100.00%
		4	0.81%
		18	0.22%
Self-developed method	13	0 (initial)	100.00%
		4	0.19%
		18	0.06%

As can be seen in the following figure, the optimum is nearly reached after a few iterations. In this steep initial part, both algorithms required approximately 2 iterations, on average, to reduce the target value by one order of magnitude.

After reaching 0.25%, the Broyden algorithm stayed stable while the new algorithm converged slowly to a lower value.

DISCUSSION OF THE RESULTS

Both algorithms converged in all of the studied cases until either the numerical accuracy limit or a limit due to an over-constrained system was reached. The number of iterations can be approximated by the formula:

(11) $\mathit{m}=1+\mathit{k}+\mathit{r}★\mathit{c}$

where:

m is the total number of iterations including the first;

k is the number of parameters (iterations of Jacobian matrix generation); and

r is a value specific to the problem; it depends mainly on the nonlinearity of the model and the interactions between the parameters.

The total number of parameters is not very significant for this value. For the load cases studied, the value was between 2 and 10, on average, until a certain limit was reached. For similar load cases, it is not expected to increase significantly; c is the desired reduction in the order of magnitude of the difference of the root sum square and the minimal possible root sum square. It is defined by the following formula:

(12) $\mathit{c}={\mathit{log}}_{10}({\mathit{R}}_0-{\mathit{R}}_{\min })-{\mathit{log}}_{10}({\mathit{R}}_{\mathit{final}}-{\mathit{R}}_{\min })$

where:

R₀ is the initial root sum square of the deviation vector (F(p_o) - t_mes);

R_final is the final root sum square of the deviation vector (F(p_final) - t_mes); and

R_min is the minimum possible root sum square of the deviation vector (F(p_final) - t_mes).

As can be seen in Fig. 6 the number of iterations needed for the presented algorithms (10 to 30) to the typical number of iterations needed for genetic (25600 (van Zijl, 2013), 843 to 33072 (Jouffroy, 2007)), adaptive particle swarm optimization (6000 (van Zijl, 2013)) or the stochastic design improvement (555, (Mareschi et al., 2005)), this approach is much faster. The number of parameters of the models used in these evaluations is comparable to the complex model used in this paper, but the number of nodes used is considerably lower. Unfortunately, the models and tools used for these methods are not accessible so that a direct comparison was not possible.

Figure 6
Comparison of the number of iterations needed for thermal model correlation (n=nodes of the model, p=parameters correlated).

For each iteration, the CPU time required for generating a new set of parameters is only a few milliseconds and, therefore, it is negligible compared to the minutes or hours required for solving the model. It can be concluded that the total CPU time is nearly proportional to the number of iterations, independent from the algorithm.

Both algorithms studied showed similar performance. Therefore, both are considered suitable for model correlation. The new algorithm showed better performance for most situations.

CONCLUSION

It has been shown that using root finding algorithms of the Broyden class for thermal mathematical model-to-test correlation is feasible and efficient. Although a direct comparison of the analyzed methods and a genetic algorithm was not possible, the number of iterations needed for the tested algorithms is 20 to 1,000 times smaller than the number of iterations reported for genetic or adaptive particle swarm algorithms. In conclusion, this algorithm significantly reduces the cost of a thermal model correlation.

The "good Broyden method" and the self-developed method delivered satisfactory results and are, therefore, both considered suitable. The self-developed method delivered slightly better results for most of the cases analyzed.

ACKNOWLEDGMENTS

The laser communication terminal (LCT) project, which this publication is based on, is carried out on behalf of the Space Administration of the German Aerospace Center (DLR e.V.) with funds from the Federal Ministry of Economics and Technology, under the reference number 50-YH1330. The author is responsible for the contents of this publication and gratefully acknowledges their support.

REFERENCES

• Broyden, C.G., 1965, "A Class of Methods for Solving Nonlinear Simultaneous Equations", AMS Mathematica of Computation, Vol. 19, No. 92, pp. 577-593. doi: 10.1090/S0025-5718-1965-0198670-6.
  » https://doi.org/10.1090/S0025-5718-1965-0198670-6
• De Palo, S., Malost, T. and Filiddani, G., 2011, "Thermal Correlation of BepiColombo MOSIF 10 Solar Constants Simulation Test", 25^th European Workshop on Thermal and ECLS Software.
• Harvatine, F.J. and DeMauro, F., 1994, "Thermal Model Correlation Using Design Sensitivity and Optimization Techniques", 24th International Conference on Environmental Systems and 5th European Symposium on Space Environmental Control Systems.
• Jouffroy, F., 2007, "Thermal model correlation using Genetic Algorithms", 21st European Workshop on Thermal and ECLS Software.
• Mareschi, V., Perotto, V. and Gorlani, M., 2005, "Thermal Test Correlation with Stochastic Technique", 35th International Conference on Environmental Systems (ICES).
• Momayez, L., Dupont, P., Popescu, B., Lottin, O. and Peerhossaini, H., 2009, "Genetic algorithm based correlations for heat transfer calculation on concave surfaces", Applied Thermal Engineering, Vol. 29, No. 17-18, pp. 3476-3481. doi: 10.1016/j.applthermaleng.2009.05.025.
  » https://doi.org/10.1016/j.applthermaleng.2009.05.025
• Roscher, M., 2006, "Genetische Methoden zur Optimierung von Satelliten Thermalmodellen bei EADS ASTRIUM", Hochschule Wismar & EADS Astrium.
• van Zijl, N., 2013, "Correlating thermal balance test results with a thermal mathematical model using evolutionary algorithms", Faculty of Aerospace Engineering, Delft University of Tecnology.
• WenLong, C., Na, L., Zhi, L., Qi, Z., AiMing, W., ZhiMin, Z. and ZongBo, H., 2011, "Application study of a correction method for a spacecraft thermal model with a Monte-Carlo hybrid algorithm", Chinese Science Bulletin, Vol. 56, No. 13, pp. 1407-1412. doi: 10.1007/s11434-010-4053-z.
  » https://doi.org/10.1007/s11434-010-4053-z

Publication Dates

• Publication in this collection
  Oct-Dec 2014

History

• Received
  26 May 2014
• Accepted
  06 Oct 2014