Induction motor control: homotopy continuation approach and power efficiency maximization


Induction motor control: homotopy continuation approach and power efficiency maximization

Alex Borisevich

June 21. 2012

Introduction




   The goal of this report is to present two main results obtained in research work during March-May 2012 in Laboratory of Electrical Drives at Reutlingen University. This work devoted to developing some advanced control algorithms for induction motors. Research is partially supported by Robert Bosch Center for Power Electronics (RBZ) and SEW-Eurodrive GmbH.
   In first chapter abstract problem of output setpoint tracking for affine non-linear system is considered. Presented approach combines state feedback linearization and homotopy numerical continuation in subspaces of phase space where feedback linearization fails. The method of numerical parameter continuation for solving systems of nonlinear equations is generalized to control affine non-linear dynamical systems. Application of proposed method demonstrated on the speed and rotor magnetic flux control in the three-phase asynchronous motor.
   In second chapter the problem of maximizing the efficiency of the induction motor under part-load condition considered. By measuring the difference between power consumption of the quadrature and direct supply channels in dq coordinate system a simple criterion for optimal power loss operation is obtained. The approach differs in that the calculation of developed criterion requires only one parameter of the motor, the the stator inductance. The value of the criterion used in the correction of the flux (magnetizing current) for its optimality tuning.
   The author is grateful for all help and constructive discussion to professor of Technical Faculty Dr.-Ing. Gernot Schullerus. Also author must thank for the technical and motivational support staff of Technical Faculty, namely Prof. Dr.-Ing. Jurgen Trost, Dr. Bernd Petereit, Dr. Daniel Fierro.

1

Chapter 1





                Switching strategy based on homotopy continuation for non-regular affine systems with application in induction motor control




            1.1    Motivation


   Let the affine nonlinear system with m inputs and m outputs in state space of dimension n is given:

m
x = f(x) + y^gi(x)ui, y = h(x),                     (1-1)
i=1
   where x е X C Rⁿ, y e Y C Rm, u e U C Rm, maps f : Rⁿ ^ Rⁿ, gi : Rⁿ ^ Rⁿ, h : Rⁿ ^ Rm are smooth vector fields f,g,h e C'■. Func tions f (x) a nd g(x) are considered as bounded on X.
   Systems of the form (1.1) are the most studied objects in the nonlinear control theory.
   There are several most famous control methods for systems of type (1.1) : feedback linearization [1, 2, 3], application of differential smoothness [4], Lyapunov functions and its generalizations [5], including a backstepping [6], also sliding control [7] and approximation of smooth dynamic systems by hybrid (switching) systems and hybrid control [8].
   All of these control techniques have different strengths and weaknesses, their development is currently an active area of research, and the applicability and practical implementation has been repeatedly confirmed in laboratory tests and in commercial hardware.
   Approach described below is based on the method of numerical parameter continuation for solving systems of nonlinear equations [9], which deals with parametrized combination of the original problem, and some very simple one with a known solution. The immediate motivation for the use of parameter continuation method in control problems is a series of papers [10, 11], in which described the application of these methods directly in the process of physical experiments.
   In this paper we consider the solution of the output zeroing problem for the system (1.1) with relative degrees rj > 1 that expands earlier obtained in [12] and [13] results for a case rj = 1. Further it is supposed that (1.1) it is free from zero-dynamics, i.e. n = 52”=₁ rj.
   The article consists of several parts. We briefly review the necessary facts about the method of parameter continuation and feedback linearization. Next, we represent the main result, an illustrative example of the method, as well as an example of controlling three-phase induction motor.


            1.2    Problem statement and motivation


   In this paper we consider the problem of nonlinear output regulation for affine nonlinear system. In particular, we will solve the problem of output regulation to constant setpoint (without loss of generality, regulation to 0).

2

Definition 1. Given the system of form (1.1). Problem of output regulation to zero (aka output zeroing) is Uie design of such state-feedback a)ntr= Zaw u(t) = u(x) application of which asymptotically drives the sy^t^m output to 0: lim, . ᵥ y(t) = 0.

    The output zeroing problem of affine nonlinear systems can be solved using mentioned above feedback linearization method. The main idea of the method consists in the transformation using a nonlinear feedback nonlinear system N : u(t) ^ y(t) to the linear one L : v(t) ^ y(t) with the same outputs y, but new inputs v. After that, the resulting linear system L can be controlled by means of linear control theory.
    Suppose that a control problem of N can be in principle solved, i.e. there is exists a satisfying input signal u*(t), which gives the output response y*(t). The essence of problems in feedback linearization comes from that the response y*(t) may not be in any way■ reproduced by system L which is obtained after linearization.
    The simplest specific example is the system x = u, y = h(x) = x(x² — 1) + 1, x(0) = 1 for which the problem of output zeroing y ^ 0 is needed to solve.
    If the system under consideration was a constant relative degree, the use of control v = —y after feedback linearization would give the output trajectory of y(t) = exp(—t), which is everywhere decreasing y(t) < 0.
    In this case, the nonlinearity y = h(x) has two limit points xj ₂ = ±1/V%, in which h'ₓ(x( ₂) = 0. Any trajectory y(t), that connects y(0) = 1 with y(T) = 0 passes sequentially through the points yj = h(3⁻¹/²) and y2 = h(—3⁻¹/²), and besides yj > yj. Hence, anу trajectory y(t) on the interval (0,ti) should decrease with time (Figure 1), on the interval (ti,t₂) increase, and in the interval again decrease. Such a trajectory is not reproducible using the feedback linearization.


Figure 1. Output trajectories of linearized system with constant relative degree and system with y = h(x) = x(x² — 1) + 1

    The behavior of the system in Figure 1 can be interpreted as follows: in the intervals (0, ti) and (t₂, T) the system can be linearized in the usual manner and presented in the form y = v. On the interval (ti, t₂) system behavior differs from the original, and the trajectory need to move in the opposite direction from the y = 0, which is the same as control of system y = —v. A similar situation arises in numerical methods for finding roots and optimization of functions with singularities, where in order to achieve optimum or find a root motion in the direction opposite to predicted by Newton’s method is needed. We can use the parameter A € [0,1] to indicate the motion direct ion. Increasing of parameter A > 0 corresponds to the movement of y(t) in the direction to the desired setpoint y = 0, and parameter decreases A < 0 in the opposite movement. The points of direction change A(t) = 0 correspond to overcoming the singularities of h(x). In fact, this idea is the basis of the approach proposed below.



            1.3  Background


  In this section we present known facts needed to understand the main result. Finally, we come to the conclusion that the numerical homotopy methods can be used not only for solving nonlinear equations, but also for control of nonlinear affine systems. Here and below we will always consider a setpoint tracking problem.

        1.3.1   Feedback linearization

Definition 2. MIMO nonlinear system has relative degree rj for оutput yj in 6 C Rⁿ if at least for one function gi is true


L Lf 1hj =0                                               (1.2)
    where Lf X = 'ff⁴ f (x) = P”=i '''^fx fi(x)     a Lie deriva tive of function X along a vector field f.
    It means that at least one input uₖ influences to output yj after rj integrations.
    Number r = 52 Г i rj is called as the total relative degree of system. If r = n and matrix

•Г1 - 1 7

A(x) =

/ LgᵢLf¹⁻¹hi(x)

L
 gm

\LS₁ Lf

hm (x)

L„ Lf
gm f

hi(x)

hm(x)

(1-3)

is full rank, then the original dynamical system (1.1) in 6 equivalent to system:

yjrj) = Lf hj + Xх LgiLf 1hj • Ui = B(x) + A(x) • u i=1

(1-4)

The nonlinear feedback

u = A(x) 1[v — B (x)]                                      (1-5)
   converts in subspace 6 original dynamical system (1.1) to linear:

y⁽rj ) = Vj                                          (1.6)
    Control of a nonlinear system (1.1) consists of two feedback loops, one of which implements a linearizing transformation (1.5), second one controls the system (1.6) by any known method of linear control theory.
    A significant drawback, which limits the applicability of the feedback linearization in practice is requirement of relative degree r constancy and full-rank of matrix A(x) in the whole phase space 6.


        1.3.2 Numerical continuation method

   Let it is necessary to solve system of the nonlinear equations

ф(£) = о

(1-7)

    where ф : Rm н Rm is vector-valued smooth nonlinear function.
    Lets Q C Rm is open set and C(Q) is set of continuous maps from its closure Q tо Rm. Functions F₀,F1 e C(Q) are homotopic (homotopy equivalent) if there exists a continuous mapping

H : Q x [0,1] н Rm

(1-8)

   that H (£, 0) = F₀(£) a nd H (£, 1) = F1(£) for all £ e Q. It can be shown [9] that the equation H(£, X) = 0 has solution (£, X) for all X e [0,1]. The objective of all numerical continuation methods is tracing of implicitly defined function H(£, X) = 0 foг X e [0,1].
   Lets H : D н Rm is C-^continuous function on an open set D C Rm⁺¹, and the Jacobian matrix DH(£, X) is full-rank rank DH(£, X) = m for all (£, X) e D. Then, for all (£, X) e D exists a unique vector т e Rm⁺¹ such as

DH(£, X) • т = 0, кт||₂ = 1, det (DH<’ X⁾"j > 0,

and mapping

Ф : D ■ Rm⁺¹, Ф : (£, X) н т

(1-9)



(1-Ю)

is locally Lipschitz on D.
Function (1.17) specifies the autonomous differential equation

                    d /А
                    - | =ф«,а), e(0) = eo, л(о) = о dt \лJ

(1-11)

    which has a unique solution ^(t), A(t) according to a theorem of solution existence for the Cauchy problem.
    It can be shown [9] that integral curve y(t) = (^(t), A(t)) reaches the solution point (£*, 1) where ф(£*) = 0 in Unite time t < to.
    In [12,13] the method of numerical homotopy continuation for the time-dependent systems of nonlinear equations ф(£, t) = 0 is given. Another variant of numerial continuation for the system of nonstationary equations described in [15].



            1.4 Main result


   In [12] pointed that the problem state control of the affine system (1.1) with relative degree of each state rj = 1 associated with the solution A" (t) : R+ ^ Rⁿ for nonstationary nonlinear equation ф(£, t) = 0, in which ^Ф = g(z), уф = f (z)- Namely, the control u(t) = d^d⁽t⁾ such that ф(£*(0),0) = x(0), brings the system (1.1) from the state x(0) tо x(T) = 0 asymptotically T ^ to. Below we give another method than that described in [12] that does not generate discontinuous control trajectories near limit points.

        1.4.1   Numerical continuation method for nonstationary system of nonlinear equations

   Consider the solution of the vector nonstationary equation ф(£, t) = 0, where ф : Rⁿ ^ Rⁿ is a smooth function. Compose parameterized simultaneously over time t and parameter Л homotopy map:


H⁽^,л,t⁾ = ⁽¹ — л⁾ • (£ — €o) + л • ф⁽£, t)

(1-12)

   where £₀ is initial approximation to the solution.
    Lets formulate and give short proofs of some assertions for the background of parameter continuation method for nonstationary system of nonlinear equations.
Assumption 1. If £(t) a nd Л(^ are solution of (1.12), t/ien rank DH^,\(^(t), X(t) ,t) = m.
    From this assumption it follows that along the trajectory (£(t), Л(£)) there is no bifurcation points in which rankDH^x(£, Л, t) < m and D\H e im/ф-H. In other words, the curve (^(t), Л(^)) is free from branches and self-intersections.
    Equation H(£, Л, t) = 0 foг Л e [0,1] defines implicitly defined function £(t), parameterized by Л(^) and satisfies the equation obtained by differentiating (1.12) by time


•                dH • (Ш • d H
HT(e(t), Л(^), t) = — £(t) + —Л(t) + — = 0                             (1.13)


   Denoting

dH dH\
1% злу

dH ~dt’

(1-14)

B =

    (1.13) can be represented as a linear matrix equation with respect to Л and Л:


At = B


(1.15)

Lemma 1. Equation (1.15) always has a solution.


Proof Because of A = DH^,\(A Л, t) the undetermined equation (1.15) has a solution if and only if rankDH^,\(£, Л, t) = m, what is based on the assumption 1. □

Lemma 2. All the solutions of (1.15) can be represented ts T = a • t + T, where T = A+B, t e ker A, a e R.


Proof It is quite obvious, since dim ker A = 1 and all the null-space of A can be parametrized by scalar variable a e R, and the solution space W inhomogeneous equation of the form (1.15) is defined as W = {A+B} ф ker A. □

   Let’s prove assertion of Lipschitz maps Ф : D н Rm+1, which will need further.

Theorem 1.     Let H : D н Rm is C¹ -continuous function on an open set D C Rm⁺2, and the Jacobian matrix A is full-rank rank A = m for all (5, A, t) e D. Then f or each (5, A, t) e D exists a unique vector т e Rm⁺¹, such that


т = a • т + т, т = A+B,

A
A • т = 0, ктЦ2 = 1, det ( T > > 0, a = const e R+,


   and map


L : D ■ Rm⁺¹, L : (5,A,t) н т

(1.16)

(1-17)

    is locally Lipschitz on D.


Proof Vector т = a • т + f defined as the sum of two components, one of which т is known (initial value problem 2.1.9 in [9]), that it is a Lipschitz function on D. Hence, it is necessary to prove that f = A+B = (DH^,\)+DHₜis Lipschitz. The uniqueness of the f follows from the uniqueness of the Moore-Penrose matrix pseudo-inversion.
    We can assume that DH is Lipschitz on D with a constant 7, which leads to existence and boundedness of the second derivatives H. Si nee rank A = m, the pseudo-inversion function A+ continuous and differentiable [11]. Hence the product A+B is Lipschitz, because its components are Lipschitz. □

    Function Ф specifies the autonomous differential equation


                             d /A ~ ,         ,   , ,       , ,                            ,
                             - 5 = L(5, A, t), 5(0) = 5o, A(0) = 0                         1.18
                             dt A
   which has a unique solution 5(t), A(t) according to a theorem of existence and uniqueness of solutions of the Cauchy problem.

Theorem 2.    The set H-¹ ({0}) is simply connected.

Proof By theorem 2.1 from [20], if for the map F : Rⁿ н Rk, k < n true that


sup{||(DF(x)DTF(x)) 1|| | x e Rⁿ} < to

(1.19)

   then F ¹ (0) connected submanifold of dimension (n — k) в Rⁿ.
   In accordance with the assumption 1, rank DH = m, therefore rank DTH = m, rank DH • DTH = m, inversion of the matrix (DH • DTH) ¹ defined and its norm is bounded. □


Theorem 3.     A necessary condition for the ex^steAce of integral curve (5(t), A(t)) for (1.18) that connects tho poonts (50,0) and (5*, 1) i s a> 0.


Proof Based on lemma 2 and theorem 1 it is obvious that the Cauchy problem of the form (1.18) has a unique solution (5(t), A(t)), that satisfies the equation (1.12) for each fixed t = const: H(5(t), A(t),t) = 0.
    To determine the sign of a consider the behavior of the curve y(t) = (5(t), A(t)) ne ar t = 0. Since


A =      + (1 — A)E ф — eY B = — ,ф
\ d5                J        di

(1.20)

    then near t = 0 equation (1.12) behaves as a stationary with B = 0. Let’s show that a = 0. If a = 0, then at time t = 0, from lemma 2 and (1.20) it follows that 5(0) = 0, A(0) = 0, thereof A(t) = 0 for each t

and condition A =1 never attainable. Since the near t = 0 equation (1.12) stationary and the solution is determined by the conventional method of parameter continuation [9], then orientation of vector т must A
be agreed with the condition ( T > > 0. He nee a > 0. □


    Now we state and give a short proof of the assumption about the behavior of the solution curve ₇(t) = (e(t),A(t))of (1.18).

Theorem 4.      There exists a number a₀ e R, that the intepraZ orrw 7(t) = (5(t), A(t)) for equation (1.18) with a > a₀ has finite length between tZie poonts (50,0) aⁿd (5*, 1)-

Proof To prove this we consider the structure of the right side (1.18). Since the constant a can be selected arbitrarily large, then the term т can be neglected and may be written:


т = a • т + т------> a • т                                 (1-21)
а^то
    By theorem 1 the function т is Lipschitz on D. then it follows automatically that т is bounded. Then one can always choose a finite a so that т = a • т + e, e C a • т.
    Since a is a finite number, then a • т as right side of (1.18) satisfies the well-known results on the finiteness of the solution trajectory (lemma 2.1.13 and theorem 2.1.14 in [6]). Hence, for an appropriate choice a> a₀ curve (£(t), A(t)) has no limit points, and is diffeomorphic to the line, i.e. has finite length between the Ao = 0 и A1 = 1. □

    The last theorem indicates that the parameter a is another one degree of freedom in designing the controller. The larger this constant, the faster the solution arrives to the A = 1, but numerical integration becomes more stiff.


        1.4.2   Homotopy continuation for nonlinear affine systems

   Let’s associate with the plant (1.1) linear dynamics system with m inputs u, n states z, m outputs n and with the same relative degries ri for outputs such as in (1.1)


z = Az + Bu, n = Cz,

d⁽rⁱ⁾
n , nn = ui.
dt r ' i

(1.22)

    Following equation is the homotopy mapping that links the outputs dynamics of the system (1.1) and (1.22):


H = (1 - A) • n + A • y = 0.

(1.23)

    By definition of the relative degree of output, each component Hi should be differentiated ri times with respect to t until it becomes an explicit function of any input u. We obtain after differentiation:


                        ri-1                                             ri-1
H ') = - X C*ni⁽rⁱ⁻fc⁾A⁽fc⁾ + (1 - A)ui + (yi - ni)A⁽rⁱ⁾ + X C*yir-kh^ k=1                                                         k=1
+A (L7 hi + X LgₖCf^hi • uk) =0, k=1
    that gives:

H⁽rⁱ⁾ = Ai,1(x, z, Л) • u + Ai,2(x, z, Л) • A⁽rⁱ⁾ + Bi(x, z, Л),

(1-24)

(1.25)

    where Л = (A, A, X,..., A⁽rⁱ⁻¹⁾), Ck are binomial coefficients.
    Considering all of the components Hi after differentiation according to the relative degrees of outputs ri it is possible to write an algebraic condition that specifies a continuous deformation of system (1.22) to (1.1).


H = A1(x, z, 21) • u + A^(x, z, 21) • A⁽rmax) + B(x, z, 21) = 0,

(1.26)

   where r.   = max{ri}, 21 = (A, A, A,..., A⁽rmax



    If B = 0, that is corresponds to the case ri = 1, f (x) = 0, any known method of parameter continuation (like predictor-corrector method [9]) can determine the trajectory u(t),A(t) such as ( u⁽J) ) = т(t), where A(t)
т(t) is tangent vector to the implicit curve H = 0, which is obtained from the linear matrix equation Ат = 0. This equation has infinitely many solutions as there are n conditions and n +1 variables. In order to uniquely identify u(t) a nd A(t) an addition al condition ||т || = 1 for length normalization of the vector т needed. In addition, to select the correct direction of the т, imposed a condition of its positive A A                                               A A A , AA.\
orientation relative to the surface H, given in the form of inequality det I T I > 0.
    Considering the general case B = 0, tangent vector т(t) needs to be augmented by term for notsta-tionarity compensating. Connected path (u(t), A(t)) for t e [0,T) starting at point (u₀,0) such as

lim A(t) = 1, H(u(t), A(t), t) = 0,


   can be generated from (1.26) as follows:


                                             / u \
^A(rmax ) J = a • T ⁺ T,
t = A+B,
.. ..      . A A\
A^ t = 0, кт||₂ = 1, det ( T > > 0,


(1-27)

(1.28)

   with the additional condition


rank A(x, z, Л) = m,                                 (1.29)
   where a > 0, a € R is a scalar constant, A+ is Moore-Penrose inverse of matrix A.
   Condition (1.29) is a standard assumption when using parameter continuation method, which corresponds to the possible existence of limit points of trajectories (u(t),A(t)) at which A1 / im A2, and the absence of bifurcation points. At the same time in some regions of phase space X x Z may be a situation where rank Ai(x, z, Л) < m, in that case, the system can not be linearized by the feedback, but the proposed method is applicable. Overcoming the bifurcation points, in which is observed A1 € im A2, also possible within the known approaches for the numerical parameter continuation (e.g., using the Lyapunov-Schmidt decomposition [9]).
   Equation (1.26) specifies the state feedback, and (1.28) specifies dynamics of the controller. According to (1.24), equations for B and A2 depend explicitly on the y, and thus implements the output feedback.


        1.4.3    Switching strategies for nonregular feedback linearization

   All plants in practice are subject to a variation of parameters. The variant of control offered in the previous section is definitely sensitive to parametric uncertainties in a control object. On the other hand, in the feedback linearization control the parameter variation in the plant can be compensated by the controller for the linearized system [16]. Let’s consider a hybrid method that combines the possibility of applying an external control loop and resistant to change in the relative degree of the system.
   System in form (1.24) with output H can be linearized by feedback if consider evalution of parameter A as observable internal dynamics. In case if we fix A⁽rmax) = const, A⁽rmax⁾ € {-1,1}, then we obtain from (1.24) following affine nonlinear system:


H = F(x, z, Л) + G(x, z, Л) • u
F (x, z, Л) = A2(x, z, Л) • Arm + B(x, z, Л)                      (1.30)
G(x, z, Л) = Ai(x, z, Л)
    It’s possible to write a nonlinear coordinate transformation v m- u transforming a nonlinear system (1.30) to linear one H⁽r⁾ = v


u = [Ai(x, z, Л)] ¹ (v — A₂(x, z, Л) • A⁽rmax⁾ — B(x, z, Л)


(1-31)

    Switching strategy can be described as follows: we start whith conventional feedback linearization of (1.30), then in areas S where det G(x, z, Л) « 0 and feedback linearization is not possible, it is necessary to switch control to the parameter continuation method implemented with (1.28).
    It should be noted that when v = 0 control (1-31) is a special case of parameter continuation strategy that applied far from limit points. If A⁽rmax⁾ = 0, then control (1-31) is generated by equations (1.28) with following time-dependent scaling:

a=

sign(A⁽r

max

_ ^⁽rmax)

A⁽rmax )

(1.32)

    Since the degeneracy of matrix Ai in feedback linearization control always leads to increasing to infinity at least one input in u, practical way to perform switching between control strategies is consider event of inputs saturation when max |ui| > uₘₐₓ. Input signals are always limited |u| < uₘₐₓ in the real world applications. This yields the following algorithm of the hybrid feedback linearization: we start

with positive sign s = 1 оf A⁽rmax⁾ = s and conventional feedback linearization, then in case of saturation of inputs we switch to homotopy continuation procedure, and after return of all input signals to its limits we flipped sign s := — s and switch back to feedback linearization. Formal procedure for calculating the control actions can be represented as follows

u := [Ai (x, z, Л)] ¹ (v — A2 (x, z, Л) • s — B(x, z, 2V)) if max |ui | > uₘₐₓ then

U u A                       ....           A A A
QXCrmax)) = T, A • T = ⁰, HT112 = 1, det J (       )) = A⁺B
i \\mnax) I

kₓ := (sign(A⁽rmax⁾) — A⁽rmax⁾)/A⁽rmax⁾

>0

k„ := (sign(max |u|)uₘax — max |и^|)/max |u^|

if kᵤ > k\ then

u := u • k\ s := sign(A⁽rmax A⁽rmax) — s

elseku > 1

u := u • kᵤ
A⁽rmax) := A⁽rmax ) • ku

end if

u = u + u, A⁽rmax⁾ = A⁽rmax⁾ + A⁽rmax⁾

end if

    Since the set of regions S forms a compact space, using of this algorithm prodice almost linearized system. The situation where x € S can be considered as a perturbation action in the output H.



            1.5 Applications


        1.5.1  One illustrative example

   Consider following abstract example of MIMO system, that changes its relative degree in the state space


3
  X1 = ui + X2
X 2 = u₂ + xi yi = x3 — xi + 1 У2 = x4 cos(2x2)

(1.33)

    with initial conditions x(0) = (1,1)T. We need to solve the problem of output zeroing y ^ 0. Let’s suppose that input signals constrained by inequality |ui| < 20.
    Differentiating the outputs, we obtain

yi = (3x2 — 1) • (ui + x³) = giiui + fi

У2 = (4x23 cos(2x2) — 2x4 sin(2x2)) • (u₂ + x³) = 522^ + /2

(1-34)

    Obviously, the system in interval x € [0,1]² can not be completely linearized by the feedback, because there are exists such x* th at gn(x*) = 0 о r g₂₂(x*) = 0.
    Let’s associate with (1.33) linear system of a form

ni = ui, 2)2 = u2

(1.35)

with initial conditions n(0) = (0,0)T
According to the equation (1.24) we obtain for H = 0 the following

Ai = fg“     °) • u + (1 — A)E, A2 = y — n, B = —
      A ⁰   522/

(1.36)

   The model in Simulink to control the system shown in figure 1. Modeling results are shown on figures 2-3.