*3.3. Combustion Fitting Approach*

The Newton-Raphson method applied for a single variable function is briefly introduced. If the initial estimate of the root is *xn*, a tangent can be extended from the point [*xn*, *f*(*xn*)]. The point where this tangent crosses the x-axis represents an improved estimate of the root. The first derivative *f* (*xn*) is equivalent to the slope of the tangent and can be derived on the basis of a geometrical interpretation (Figure 6) and an iteration scheme is constructed:

$$f'(\mathbf{x}\_n) = \frac{f(\mathbf{x}\_n) - 0}{\mathbf{x}\_n - \mathbf{x}\_{n+1}} \tag{2}$$

Then

$$\mathbf{x}\_{n+1} = \mathbf{x}\_n - \frac{f(\mathbf{x}\_n)}{f'(\mathbf{x}\_n)} \tag{3}$$

While |*xn+*1 − *xn+*1| ≤ *tolerance*, the *xn+*1 can be considered to be the root of function *f*(*x*).

**Figure 6.** Graphical description of Newton-Raphson method.

Furthermore, the Newton-Raphson method can be used to acquire the solution for a multi-variable function **f(x)**. If the zero was considered to occur at **x** = **<sup>x</sup>\***, where **x\*** is a vector, then the Taylor series for multi-variable function is applied. If **xn** is assumed to be the current estimation of the functions and **xn+1** = **xn** + **h**, due to the demand that **f**(**x\***) = **0**, the Taylor series is acquired:

$$\mathbf{f(x)}|\_{\mathbf{x}=\mathbf{x\_{n+1}}} = \mathbf{f(x\_n+h)} = \mathbf{f(x\_n)} + \nabla \mathbf{f}|\_{\mathbf{x\_n}} \cdot \mathbf{h} + \mathbf{o}\left(|\mathbf{h}|^2\right) = \mathbf{0} \tag{4}$$

οIf **o**-|**h**|<sup>2</sup> is assumed small enough to neglect, Equation (4) becomes:

$$\mathbf{f(x\_n)} = -\nabla \mathbf{f}|\_{\mathbf{x\_n}} \cdot \mathbf{h} \tag{5}$$

When − <sup>∇</sup>**<sup>f</sup>**|**xg**is replaced by matrix **A**:

$$\mathbf{h} = -\mathbf{A}^{-1} \cdot \mathbf{f}(\mathbf{x\_n}) \tag{6}$$

$$\mathbf{A} = \begin{bmatrix} \frac{\partial f\_1}{\partial x\_1} & \frac{\partial f\_1}{\partial x\_2} & \cdots & \frac{\partial f\_1}{\partial x\_n} \\ \frac{\partial f\_2}{\partial x\_1} & \frac{\partial f\_2}{\partial x\_2} & \cdots & \frac{\partial f\_2}{\partial x\_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f\_n}{\partial x\_1} & \frac{\partial f\_n}{\partial x\_2} & \cdots & \frac{\partial f\_n}{\partial x\_n} \end{bmatrix} \text{and } \mathbf{f}(\mathbf{x}) = \begin{bmatrix} f\_1(\mathbf{x}\_1, \mathbf{x}\_2, \cdots, \mathbf{x}\_n) \\ f\_2(\mathbf{x}\_1, \mathbf{x}\_2, \cdots, \mathbf{x}\_n) \\ \vdots \\ f\_n(\mathbf{x}\_1, \mathbf{x}\_2, \cdots, \mathbf{x}\_n) \end{bmatrix} \tag{7}$$

Finally, the next estimate of **xn+1** is obtained and used as the initial value of the next iteration and the iteration is terminated when the stop criteria are met.

### *3.4. Implementation in the Combustion Fitting Based on Seiliger Process*

The Newton-Raphson multi-variable root finding method is applied to find the solution for the fitting of the (real) engine equivalence values by means of varying Seiliger process parameters. The Seiliger process is an engine combustion model to divide the engine working cycle into finite stages based on different Seiliger parameter numbers choosing. The case with four equivalence criteria as functions of four variables (Seiliger parameters) is taken as an example.

Seiliger parameters *a*, *b*, *c* and *nexp* are the variables in these equations and the differences between the calculated result from the equivalence criteria functions of the Seiliger process (*pmax*, *Tmax*, *qin*,*Seiliger* and *wi*,*Seiliger*) and the measured engine cycle(*p*3, *T*4, *qin*,*measured* and *wi*,*measured*) are the functions for which the zero has to be found. The other parameters in Equations (8)–(11) are considered to be constant but changed with engine working conditions. Therefore, the simulation model is important to determine them at different iteration steps rather than analytical method.

In summary, the Equations (8)–(11) are expressed by standard equation form

$$f\_1(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3, \mathbf{x}\_4) = p\_3 - p\_{\max} \tag{8}$$

$$f\_2(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3, \mathbf{x}\_4) = T\_4 - T\_{\max} \tag{9}$$

$$f\_3(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3, \mathbf{x}\_4) = q\_{\text{in\\_Sailiger}} - q\_{\text{in\\_measured}} \tag{10}$$

$$f\_4(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3, \mathbf{x}\_4) = w\_{i, \text{Senliger}} - w\_{i, \text{measured}} \tag{11}$$

The five functions are set to a column vector as **func**:

$$\mathbf{\texttt{runc}} = \begin{bmatrix} f\_1 \ \ f\_2 \ \ f\_3 \ \ f\_4 \end{bmatrix} \tag{12}$$

A matrix **PDE** then is set to solve the derivatives:

$$\mathbf{PDE} = \begin{bmatrix} \frac{\partial f\_1}{\partial x\_1} & \cdots & \frac{\partial f\_1}{\partial x\_4} \\ \vdots & \ddots & \vdots \\ \frac{\partial f\_4}{\partial x\_1} & \cdots & \frac{\partial f\_4}{\partial x\_4} \end{bmatrix} \tag{13}$$

$$\mathbf{Vect} = -\mathbf{PDE}/\mathbf{func}\tag{14}$$

**Vect** is the increment or decrement of the root that is obtained from the previous step:

$$\mathbf{x}\_{i} = \mathbf{x}\_{i} + \mathbf{Vect}\ (i) \quad \text{1, 2, 3, 4} \tag{15}$$

The partial derivatives of the elements in matrix **PDE** are obtained with numerical analysis rather than with an analytical method. The latter is impossible because some parameters in equivalence criteria functions are interdependent and a numerical method is used to avoid this problem.

Variable *xi* is changed one percent and then the partial derivatives are obtained relying on the local linearity.

$$
\Delta \mathbf{x}\_i = 0.01 \mathbf{x}\_i \quad i = 1, 2, 3, 4, 5 \tag{16}
$$

$$\frac{\partial f\_i}{\partial x\_j} = \frac{f(\mathbf{x}\_{\text{'}} + \Lambda \mathbf{x}\_{\text{'}}, \dots) - f\left(\mathbf{x}\_{\text{'}}, \dots\right)}{\Lambda \mathbf{x}\_j} \quad \mathbf{i} = 1, 2, 3, 4, 5; \quad \mathbf{j} = 1, 2, 3, 4, 5 \tag{17}$$
