# Solving a Class of Nonlinear Inverse Problems Using a Feedback Control Approach

Authors: Mahmoud Tarokh
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Inverse problems have applications in many branches of science and engineering. In this paper we propose a new approach to solving inverse problems which is based on using concepts from feedback control systems to determine the inverse of highly nonlinear, discontinuous, and ill-conditioned input-output relationships. The method uses elements from least squares solutions that are formed within a control loop. The stability and convergence of the inverse solution are established. Several examples demonstrate the applicability of the proposed method.


---

## Body

## 1. Introduction

Inverse problems form an interdisciplinary filed encompassing science and engineering. Applications of inverse problems are diverse and include robotics, optics, nondestructive detection, geophysics, imaging, acoustics, and civil and mechanical engineering, just to list a few. Briefly and roughly speaking, “forward” problems estimate the effect (output) from the cause (input). In contrast, “inverse” problems require estimating the cause or the parameters from the effect. For example, in robot manipulators the forward kinematics equations are readily available and relate the robot manipulator joints angles to the end effector (hand) position and orientation. However the fundamental problem in robotics is the inverse kinematics; that is, given hand position and orientation, determine the joint angles. The mathematical formulation of inverse problems leads to models that are typically ill-posed, meaning that a solution may not exist, the solution may not be unique, or the solution is sensitive to small changes in data causing instability. In the robot manipulator example, this situation often happens when during the robot motion the joint angles form a configuration in which the so-called manipulator Jacobian matrix is singular.Inverse problems can be categorized into two groups: linear and nonlinear. Linear inverse problems are usually formulated asy=Ax where y is vector of (possibly noisy) data and the input vector x is to be extracted, given that the matrix A can be ill-conditioned such that its inverse either does not exist or is meaningless and does not reflect the actual physical problem. In other to overcome these difficulties, one has to use regularization which replaces an ill-posed problem by a neighbouring well-posed one. There are two main related regularization approaches, namely, Tikhonov and truncated singular value decomposition. In the Tikhonov approach a small constant is included in the pseudo-inverse computation of the matrix A to limit to magnitude of x. In the truncated SVD method, tiny singular values of A are removed to prevent the pseudo-inverse of A to become excessively large. Variations and improvement of these two basic approaches have been proposed (e.g., [1–3]).Many real-world inverse problems are nonlinear and unlike the linear ones have not been fully explored due to the complexity of the problem. Many of the existing techniques are based on least squares formulation of the ill-posed nonlinear problem of the formy=f(x) [4–8]. Most damped least squares solutions are based on Levenberg–Marquardt iterative solution which is further explored in [9]. A more recent work describes a two-stage method combining Levenberg–Marquardt regularization of the linearized problems at each stage [10]. When the solution is assumed to have a sparse expansion, an iterative method has been proposed to solve the nonlinear inverse problem [11]. This method shows good results for colour image inpainting applications.A recent book describes the application of inverse problems to imaging [12]. A survey of nonlinear inverse problems can be found in [13] and references therein. Neural networks have also played a major role in the solution of inverse problems (see, e.g., [14] and references therein). This is due to the fact that neural networks are considered universal approximators. A method using function decomposition and approximation is proposed by the author that is applicable to a wide-class of nonlinear inverse problems [15].Many approaches to inverse problems propose solutions for specific and often simple cases and thus are not generally applicable to situations and applications other than those that they are developed for. In this paper we present a new approach to the inverse problems that have complex structures with highly nonlinear and discontinuous behaviours as well as being ill-conditioned. We employ concepts from feedback control systems and stability theory to prove the convergence and accuracy of the solutions. Examples are provided to demonstrate the concept.


## 2. Statement of the Problem

Consider the input vectorx of dimension n×1 and the corresponding output vector y of dimension m×1 and the forward relationship between x and y of the form (1)y=fx.It is noted that (1) also covers cases where the input-output relationship is expressed in the integral form(2)y=∫abgx,udu.We assume that the forward relationship of the form (1) or (2) is known. The inverse problem is given y; find x. In most applications there is no analytical solution to the inverse problem; that is, x cannot be expressed in analytical form in terms of y.There are a number of challenges to the inverse problems. Noise can corrupt the data, the forward and/or the inverse relationship can be discontinuous, and the problem can be ill-conditioned. Furthermore, in some applications the forward relationship can be very complex and the inverse can have multiple solutions for a particular value ofy. These are especially the case for many practical problems such as robotic inverse kinematics problem. In such an application even the forward equations consist of highly nonlinear trigonometric functions with sever discontinuities due to joint limits and other physical constraints, as well as varying number of the inverse solutions.We propose to solve the above inverse problem using a feedback control approach and a generalized pseudo-inverse method, as well as optimizing a performance index in some cases. Suppose that the desired output vector isyd and we want to find an input xd such that yd=f(xd). Suppose further that, for an initial value x0, for example, randomly chosen, the corresponding output is obtained from (1) as y0. The problem is then to transfer from y0 to yd. We propose to do this transition through a smooth reference trajectory parameterized by time yr(t). Among other possible functions such as quintic, we propose a cycloid function defined by (3)yrt=y0+tτ-12πsin⁡2πτtyd-y0t≤τydt>τ,where τ is the time to move from the initial value y0 to the desired value yd. It is noted that the first, second, and higher order time derivatives of (3) are all continuous, smooth, and finite which is an important characteristic of this particular trajectory. Higher values of τ correspond to slowly varying trajectory and have a favorable effect on the stability of the system, as will be seen later. The cycloid function has a simple form and is parameterized by only three quantities, namely, y0, yd, and τ, each of which has a meaningful interpretation. On the other hand, other smooth functions are described by more parameters whose roles in the overall form of the function are not easily recognizable. For example, the quintic function requires six parameters, and it is not easy to determine these six parameters to achieve the desired transition time, final value, and maximum value of its derivative, all of which are needed later in Sections 3 and 4. Furthermore, a simple step function is undesirable since its first and higher order derivatives are infinite, resulting in instability and large steady-state error for this application, as will be seen in Sections 3 and 4. In addition, a step function has no parameters to adjust the transition time between initial and final values.Since the problem is nonlinear we use an incremental approach by taking the time derivative of (1) to obtain (4)y˙=∂f∂xx˙=Qxx˙,where Qx is the m×n Jacobian matrix. The partial derivative is computed numerically. When the trajectory passes through the points where the function fx is discontinuous, some or all elements of Qx become infinity at certain values of x. We assume that the desired xd is not at the points of discontinuity. Consequently at the points where the derivative is infinite, we set the corresponding values of the elements of Qx to a large but finite value to prevent numerical difficulties when the inverse of Qx is found. Since we seek values of xd for finite y=yd, therefore high values of the derivative can only happen during the transition, that is, at times t<τ. As a result, even though x can be very different from the desired xd during the transient, it will become arbitrarily close to xd at steady-state, that is, t≥τ, as will be shown in Section 3.Equation (4) is a differential equation relating the incremental change in the input to the incremental change in the output. The generalized inverse solution to (4) is [16] (5)x˙=Q#xy˙+βIn-Q#xQxgx,where Q#(x) is the n×m pseudo-inverse of Qx, In is the n×n identity matrix, g(x) is an arbitrary free n×1 vector, and β is an arbitrary positive scalar. The free vector g(x) and scalar β will be used for optimization as will be discussed in Section 4. The second term on the right hand side of (5) is orthogonal to the first term. It is possible to substitute for y˙t=y˙rt from (3) into (5) to obtain x˙t=x˙rt and then integrate x˙rt to get xrt whose steady-state value xd inverse corresponds to the desiredyd. However, integration of (5) is not a viable solution due to (i) integration drift and numerical difficulties and (ii) ill-conditioned matrix Qx. This means that the matrix Qx can have singular values close to zero and as a result the inverse or pseudo-inverse of the matrix, that is, Q#(x), can produce solutions that are grossly in error, especially when noise is present. To overcome (i) we propose to form an error defined by(6)et=yrt-ytand apply a feedback signal Ke(t), where K is the m×m gain matrix. To remedy (ii) we decompose Q(x) into two matrices as follows:(7)Qx=Q1x+Q2x,where Q1x is well-conditioned with a desired conditioning number and Q2 is a residual ill-conditioned matrix that can be found in a number of different ways, two of which will be described shortly. Thus, (5) is modified to (8)x˙=Q1#xyr˙t+Ke+βIn-Q1#xQ1xgx.We will prove that the steady-state solution to dynamic equation (8) will be inverse solution xd for a desired yd when a suitable feedback matrix K is determined. The feedback block diagram representing (2) and (8) is shown in Figure 1.Figure 1
Block diagram of the feedback system for the inverse problem.It is noted that since in the proposed control schemey which has the dimension m is the quantity to be controlled, and the number of adjustable variables x is n, then for controllability we must have n>m. When n>m, the n-m additional degrees of freedom can be used for an optimization task as will be seen in Section 4. When n=m, Q1#Q1=In and the second term in (8) vanishes and it becomes(9)x˙=Q1#y˙r+Ke,where the arguments of the quantities in (9) have been dropped for convenience. The inverse problem in the above control system setting may be stated as follows. Design the controller K such that the control system is stable and the error (6) tends asymptotically to a sufficiently small value ϵ; that is,(10)limt→∞⁡et=limt→∞⁡yrt-fxt=yd-fxd=ϵgiving yd≈fxd where the bound on ϵ will be established in Section 3.We now discuss how the decomposition in (7) can be achieved. The first method is based on assigning the small singular values of Q to Q2 so that the well-conditioned Q1 is a truncated version of Q. Let the singular value decomposition (SVD) of the m×n matrix Q be expressed as (11)Q=USVT=∑i=1nσiuiviT,where U=u1,u2,…,um is a unitary m×m matrix and u1,u2,…,um are m×1 vectors and have the property UTU=UUT=Im. Similarly V=v1,v2,…,vn is a unitary n×n matrix, with the same property as U. The singular value matrix is a diagonal matrix σ=diag(σ1,σ2,…,σn), where σ1>σ1>⋯>σn≥0 are the singular values. When the matrix is ill-conditioned one or more singular values are very small or even zero. This situation causes the inverse of the matrix Q to produce erroneous results. One common method to deal with this problem is to decompose Q and separate the terms involving small singular values in (10) as [3] (12)Q=∑i=1pσiuiviT+∑i=p+1nσiuiviT≡Q1+Q2,where the second summation on the right hand side of (12) contains small singular values. In the truncated singular value decomposition method, the matrix Q is replaced by Q1which is not ill-conditioned and whose pseudo-inverse is (13)Q1#=∑i=1p1σiviuiT.The amount of truncation, that is, the value of, can be determined by specifying a desired conditioning number Cd for Q1 which is the ratio of the largest singular value to the smallest singular value of Q1; that is, (14)Cd=σ1σp.Thus, given Cd and the known maximum singular value σ1, the minimum singular σp and thus the truncation index p are found.The second method to obtainQ1 is via Tikhonov method; that is, (15)Q1#=QTQQT+λ2Im-1,where λ2 is a small scalar. It can be shown that [3] (16)Q1#=∑i=1nσiσi2+λ2viuiT.Thus, (17)Q1=∑i=1nσi2+λ2σiuiviT=∑i=1nσiuiviT+∑i=1nλ2σiuiviT.Since we defined Q in (7) as Q=Q1+Q2, and the first term on the right hand side of (17) is Q, we have for the Tikhonov method (18)Q2=-∑i=1nλ2σiuiviT.It is noted that the conditioning number of Q1 obtained from the Tikhonov method is (19)Cd=σn2+λ2σ12+λ2σ1σn.Since σ1 and σn are known, λ2 can be found for a desired conditioning number Cd. For the case of single input single output systems, there will be only one singular value which cannot be truncated and in this case only Tikhonov method will be used. It is important to note that Q1obtained in (13) by the truncated singular values and Q1 found from (16) using the Tikhonov method are not necessarily the same. However, using the above procedures, both will be well-conditioned with a desired conditioning number.Having obtainedQ1 by either method, we premultiply both sides of (9) by Q1, and noting that Q1Q1#=I, we obtain(20)Q1x˙=y˙r+Ke.Substituting Q1 from (7) into (20) we obtain (21)Q-Q2x˙=y˙r+Ke.SinceQx˙=y˙ by (4), we have(22)y˙-y˙r=Ke+Q2x˙.We substitute for x˙ from (9) into (22) and noting from (6) that e˙=y˙r-y˙, we obtain the error dynamics as(23)e˙=-Im+Q2Q1#Ke-Q2Q1#y˙r.We set the controller gain K to(24)K=Im+Q2Q1#-1k,where k>0 is a positive scalar gain. Now (23) becomes(25)e˙=-ke-Q2Q1#y˙r.The above equation represents the dynamics of the error. We will discuss the effect of the gain k on the performance of the system in Section 3.To summarize this section, we have formulated the nonlinear inverse problem using incremental change and its associated Jacobian matrix. When the problem is ill-conditioned the Jacobian matrix will have singular values close to zero. The Jacobian matrix is expressed as the sum of two matrices consisting of an ill-conditioned matrix and a well-conditioned matrix. Employing a feedback control approach that uses only the well-conditioned part of the Jacobian matrix, the problem has been transformed into a dynamic system that describes the error dynamics, that is, the difference between the reference and actual output. In the next section, we prove that the error dynamic is stable and that the error can be made sufficiently small, thus providing an accurate inverse solution.


## 3. Stability and Convergence of Solutions

In this section we show that the dynamic equation (25) is asymptotically stable, meaning that the error goes to small values allowing the inverse xd corresponding to yd be determined precisely in the steady sate at the output of the system of Figure 1.In order to establish the stability of the nonlinear error dynamics (25) and its convergence, we consider the Lyapunov function candidate(26)V=12eTewhose derivative using (25) is(27)V˙=eTe˙=-keTe+eTQ2Q1#y˙r≤-ke2+Q2Q1#y˙rmaxe=-2kV+aτ2V,where · denotes the 2-norm and (y˙r)max=yd-y0/τ is the maximum value of y˙r which is obtained from (3), where Δy=yd-y0 is the difference between initial and desired values of the output. The 2-norm of a matrix is equal to its maximum singular value, and thus for the truncated SVD, with Q1 and Q2 given by (11) and (12), we have Q1#=1/σp, Q2=σp+1 and thus a in (27) is(28)a=Q2Q1#yd-y0=σp+1σpyd-y0.Note that σp+1/σp<1 and therefore a<yd-y0≡Δy. For the Tikhonov method with Q1# and Q2 given by (16) and (18), respectively, we obtain Q1#=σ1/σ12+λ2<1 and Q2=λ2/σ1≤1 if we set λ2≤σ1. As a result a=σ1/σ12+λ2λ2/σ1yd-y0<Δy. Thus, in both cases of truncated singular values and Tikhonov methods, a<Δy and (27) becomes(29)V˙≤-2kV+Δyτ2V.Note from (26) that V=1/2e2≡1/2W2, where W=e. Substituting W for V in (29), we obtain W˙≤-kW+Δy/τ whose solution is W≤W0-Δy/kτε-kt+Δy/kτ, where W0=e0 is the initial value. Thus,(30)e≤e0-Δykτε-kt+Δykτ.It is seen that the magnitude of the error decreases exponentially with time and its value is finite at all times. The steady-state value of the error, that is, limt→∞⁡e(t), is (31)e≤Δykτ.It is noted that since Δy=yd-y0, the closeness of the initial output values y0 to the desired value yd reduces the steady-state error. It is therefore advisable to assign several initial values x0 into (1) and pick the resulting y0 that is closest to yd. Regardless of this, the steady-state error can be made as small as desired by increasing the gain k and the trajectory time τ. Observe that large value of the gain k is applied to the system (Figure 1) that contains the well-conditioned inverse of Q1 and does not therefore create a problem. Furthermore, higher values of the gain k will also make the decay of the error faster as seen from (30). The steady-state value of x, that is, xd, which is obtained in the output of the control system of Figure 1, is the required inverse solution. Furthermore, when Q is regular (not ill-conditioned), then Q2=0. Consequently (25) reduces to e˙=-ke whose solution is (32)e≤e0ε-ktand thus the steady-state value of the errors is ess=0.When noiseξ is present, (4) becomes as y˙+ξ=Qxx˙. Assuming a bound ϵ exists on the noise such that ξ≤ϵ, then the truncation index p in (12) is chosen by the “discrepancy principle” such that ξ≤ϵη, where η>1. The details of discrepancy principle are provided in [3, 17]. In this case the truncation amount increases (or equivalently p decreases) with ϵ. In other words, the higher is the value of the noise norm ϵ, the fewer singular values will be included in Q1. A similar procedure is used to determine λ2 in the Tikhonov method. In this case, λ2 increases with an increase in ϵ [17]. Furthermore, it is well known that, in a closed loop system, in general feedback reduces that effect of noise on the output of the system, that is, x in Figure 1. To demonstrate this, we consider the above noise model. The error dynamics (25) with noise changes to e˙=-ke-Q2Q1#y˙r+ξ and considering the above bound on the magnitude of noise and following the procedure leading to (31), it can readily be shown that (31) is modified to ess≤Δy+ϵ/kτ. Thus, the effect of the noise on the error can be reduced by increasing the gain and trajectory final time.A note is appropriate in the situation where the inverse problem has multiple solutions for a given value ofyd. In such a case, the particular solution that the above method finds depends on the initial values x0 or y0. If the number of solutions is not too many, various solutions can generally be found by assigning random initial values within the range of interest and applying the method to obtain these various solutions. This aspect will be explored in Section 5.In summary, in this section using the Lyapunov stability approach, we have proven that the dynamics of the error, that is, the difference between the desired outputyd and the actual output y, is stable and that the steady-state error is sufficiently small, yielding an accurate inverse solution xd corresponding to the desired yd.


## 4. Optimization

For the case whenn>m and rankQ=m,there are (n-m) degrees of freedom that can be used for optimization of a performance index ψ(x). For example, if the magnitude of x is to be minimized we set ψx=x. Alternatively one might be interested to maximize signal to noise ratio.Suppose that the functionψ(x) is to be minimized, in which case we set the free vector g in (5) as(33)g=-∂ψx∂x.In order to show that ψx is minimized, we take its time derivative as(34)ψ˙=∂ψx∂xTx˙=-gTQ1#y˙r+Ke+βIn-Q1#Q1g,where x˙ has been substituted from (5). Furthermore, (34) can be manipulated to obtain(35)ψ˙≤-βgTIn-Q1#Q1g+gTQ1#y˙rmax+Kemax.The matrix In-Q1#Q1 is both symmetric and idempotent and therefore In-Q1#Q1=In-Q1#Q1In-Q1#Q1T. As a result the first term in (35) is negative definite and can be made as large as desired by the free parameter β. Since y˙rmax=1/τyd-y0=1/τΔy, this term can be made small by choosing a high trajectory time constant τ. Finally, the term involving the error is finite and small, as shown before. Considering these, the right hand side of (35) becomes negative implying that performance index ψ(x) decreases monotonically and the minimization is achieved. For maximization of the performance index, the free vector g is chosen as g=∂ψ(x)/∂x.Finally we observe that all developments in Sections2–4 are applicable to a linear relationship of the form y=Ax, in which case Q=A will be a constant matrix.


## 5. Examples

In this section we present two examples to demonstrate the features of the proposed method for solving challenging inverse problems.


### 5.1. Example 1

We consider the Gamma function which is described by(36)y=∫0∞zx-1εzdz.As shown in Figure 2, Gamma function is highly nonlinear and discontinuous and has multiple solutions for a given y and therefore provides a challenge as an inverse problem. This problem is ill-conditioned since (i) for a given y, the corresponding x is not unique, for example, in the range of  -5≤x≤4, it can have up to six solutions, and (ii) the solution does not change continuously with the initial condition since a small change in the initial condition can change the solution substantially due to the jumping from one function branch to another as is evident from Figure 2.Figure 2
The Gamma function showing a highly nonlinear and discontinuous function.As seen from (31) the steady sate error is inversely proportional to the product of the gain k and the trajectory time τ and is proportional to the difference between the desired and initial outputs: Δy=yd-y0. Since Δy is known, we can determine k and τ for a desired bound on error. It is also noted that the rate by which the errors decay with time is dependent on the gain, as seen by (30) and (32). Higher values of the gain produce faster response but can also produce overshoot and possible instability in discontinuous systems. These observations provide a general guideline for the selection of the gain and time constant. We consider the following three cases for all of which we choose the trajectory time in (3) as τ=300 and feedback gain k=0.2. This gives =60, and thus the error per unit of Δy becomes less than 1/60 = 0.016 when the system is ill-conditioned. Note that due to the very conservative bounds that were established in obtaining (31), the actual errors will be much smaller than those established in (31), as will be seen in the following cases. The set of values for k and τ is not unique, and other values satisfying the desired steady-state errors and decay time can be chosen. The errors tend to zero when the system is not ill-conditioned, as seen from (32). The regularization parameter in (15) is chosen as λ2=0.01 and τ=300 when the trajectory is near or enters the singularity (discontinuity).Case 1.
Let the desired output beyd=4.5, and suppose the initial value is chosen as x0=3 which corresponds to y0=2. This point is located in the right most branch at point A in Figure 2. The trajectories of the desired output yr(t), actual output y(t), the error e(t), and the acquired x(t) are shown in Figure 3. It is seen that the error decreases to practically zero. The final desired value of x(t) is xd=3.7646. This solution is located to the right of point A which is the closest to the initial point.Figure 3
Case1: the reference trajectory yr(t) (a), actual trajectory y(t) (b), error e(t) (c), and inverse solution x(t) (d). The maximum error is emax=6×10-4 and steady-state error is ess = 0.
(a)
(b)
(c)
(d)Case 2.
Let us now choose the initial value ofx0=0.8 with the corresponding x0=1.1642 and the same desired yd=4.5 as in Case 1. The trajectories (not shown) are similarly smooth and xd=0.2038. It is noted that now the solution is to the left of the initial point A which is the closest point to the initial value, that is, minimum change in x.Case 3.
Now suppose we use the same initial values as in Case1: that is, x0=3, y0=2, which is point A in Figure 2. However, we now change the desired value to yd=-3.5. As seen from Figure 2 this requires the passage through one or more discontinuities/singularities. The resulting trajectories are now shown in Figure 4. It is seen that, at about the time t≈100, the value of x is x≈1 which is close to the minimum value of y which is y≈1, that is, at point B in Figure 2. As a result the system oscillates to the right and left of the minimum point B to reduce the error. However, values of y around B are all positives, and achieving yd=-3.5 is impossible. As a result, the error increases until at about the time t≈140 (Figure 4) and the feedback control causes a jump to point C in Figure 2. The error now decreases until x reaches xd=-2.949 which corresponds to the desired yd=-3.5 with negligible steady-state error. The amplitude of jump depends on the gain k. Lower values of gain will make the jump to a point right of C, and higher gains make the jump to the left of point C, both of which are close to the desired yd=-3.5. This was verified and in each case an almost exact solution was found.
The above case studies and many more carried out and not reported here demonstrate the effectiveness of the proposed method in dealing with ill-condition cases, including severe nonlinearities and discontinuities.Figure 4
Case3: the reference trajectory yr(t), (a), actual trajectory y(t), (b), error e(t) (c), and inverse solution x(t) (d).
(a)
(b)
(c)
(d)


### 5.2. Example 2

We consider the complex two-input, two-output system described by(37)y1=sin⁡xcos⁡xy2=x12+2x22-x12x2;-3<x1,x2<3.The plots of the above functions are shown in Figure 5. For a desired set of (y1d,y2d), there can be multiple inverse solutions, and arriving at a particular solution depends on the initial values.Figure 5
Three-dimensional views of the functions (a)y1=sin⁡xcos⁡x and (b) y2=x12+2x22-x12x2.
(a)
(b)Case 1.
Lety1d,y2d=(0.174,1.22), x10,x20=(0.2,-0.5), y10,y20=(-0.4688,0.5600), and τ=20,k=0.5, and use the truncated singular value method to remove one singular value when it becomes below 0.01. The responses are shown in Figure 6, and the steady-state value of the output is obtained exactly equal to the desired values, and the corresponding input vector is x1d,x2d=(1.2,0.5).
We have computed the inverse desiredx1d and x2d for all values of y1d in the range of −20 to 60 and y2d from −1 to 1. The plots of inverses x1d and x2d against y1d and y2d are shown in Figure 7.Figure 6
Case1, two-input, two-output system: references trajectoriesy1r(t), y2r(t) (a), actual trajectories y1(t), y2(t) (b), errors e1(t), e2(t) (c), and inverse solution x1(t), x2(t) (d).
(a)
(b)
(c)
(d)Figure 7
Three-dimensional views of the inverse functions. (a)x1d as a function of y1d and y2d; (b) x2d as a function of y1d and y2d.
(a)
(b)Case 2.
The purpose of this case is to demonstrate the optimization discussed in Section4. Consider the two-input, single-output system of second equation in (37); that is, (38)y=x12+2x22-x12x2.Let the desired value be yd=1.22 and the initial values of x10,x20=(0.2,-0.5), as before. The corresponding initial output is y0=0.56. The optimization criterion is chosen as the minimization of ψ(x)=1/2(x12+x22) while achieving the desired output. Thus, gx=-∂ψx/∂x=-x1x2, and the parameters are chosen as τ=10,k=0.5, and β=1. The trajectories are shown in Figure 8 and indicate a smooth transition to the final values of x1d,x2d=(0.272,-0.7387) with the ψx=0.619. It is noted that with the same initial values but without optimization ψx=10.23. Thus, the optimization has minimized the change in x1,x2 while achieving the desired value.Figure 8
Case2, two-input one-output system with optimization: the reference trajectory yr (a), actual trajectory y(t) (b), error e(t) (c), and inverse solutions x1(t), x2t (d).
(a)
(b)
(c)
(d)


## 5.1. Example 1

We consider the Gamma function which is described by(36)y=∫0∞zx-1εzdz.As shown in Figure 2, Gamma function is highly nonlinear and discontinuous and has multiple solutions for a given y and therefore provides a challenge as an inverse problem. This problem is ill-conditioned since (i) for a given y, the corresponding x is not unique, for example, in the range of  -5≤x≤4, it can have up to six solutions, and (ii) the solution does not change continuously with the initial condition since a small change in the initial condition can change the solution substantially due to the jumping from one function branch to another as is evident from Figure 2.Figure 2
The Gamma function showing a highly nonlinear and discontinuous function.As seen from (31) the steady sate error is inversely proportional to the product of the gain k and the trajectory time τ and is proportional to the difference between the desired and initial outputs: Δy=yd-y0. Since Δy is known, we can determine k and τ for a desired bound on error. It is also noted that the rate by which the errors decay with time is dependent on the gain, as seen by (30) and (32). Higher values of the gain produce faster response but can also produce overshoot and possible instability in discontinuous systems. These observations provide a general guideline for the selection of the gain and time constant. We consider the following three cases for all of which we choose the trajectory time in (3) as τ=300 and feedback gain k=0.2. This gives =60, and thus the error per unit of Δy becomes less than 1/60 = 0.016 when the system is ill-conditioned. Note that due to the very conservative bounds that were established in obtaining (31), the actual errors will be much smaller than those established in (31), as will be seen in the following cases. The set of values for k and τ is not unique, and other values satisfying the desired steady-state errors and decay time can be chosen. The errors tend to zero when the system is not ill-conditioned, as seen from (32). The regularization parameter in (15) is chosen as λ2=0.01 and τ=300 when the trajectory is near or enters the singularity (discontinuity).Case 1.
Let the desired output beyd=4.5, and suppose the initial value is chosen as x0=3 which corresponds to y0=2. This point is located in the right most branch at point A in Figure 2. The trajectories of the desired output yr(t), actual output y(t), the error e(t), and the acquired x(t) are shown in Figure 3. It is seen that the error decreases to practically zero. The final desired value of x(t) is xd=3.7646. This solution is located to the right of point A which is the closest to the initial point.Figure 3
Case1: the reference trajectory yr(t) (a), actual trajectory y(t) (b), error e(t) (c), and inverse solution x(t) (d). The maximum error is emax=6×10-4 and steady-state error is ess = 0.
(a)
(b)
(c)
(d)Case 2.
Let us now choose the initial value ofx0=0.8 with the corresponding x0=1.1642 and the same desired yd=4.5 as in Case 1. The trajectories (not shown) are similarly smooth and xd=0.2038. It is noted that now the solution is to the left of the initial point A which is the closest point to the initial value, that is, minimum change in x.Case 3.
Now suppose we use the same initial values as in Case1: that is, x0=3, y0=2, which is point A in Figure 2. However, we now change the desired value to yd=-3.5. As seen from Figure 2 this requires the passage through one or more discontinuities/singularities. The resulting trajectories are now shown in Figure 4. It is seen that, at about the time t≈100, the value of x is x≈1 which is close to the minimum value of y which is y≈1, that is, at point B in Figure 2. As a result the system oscillates to the right and left of the minimum point B to reduce the error. However, values of y around B are all positives, and achieving yd=-3.5 is impossible. As a result, the error increases until at about the time t≈140 (Figure 4) and the feedback control causes a jump to point C in Figure 2. The error now decreases until x reaches xd=-2.949 which corresponds to the desired yd=-3.5 with negligible steady-state error. The amplitude of jump depends on the gain k. Lower values of gain will make the jump to a point right of C, and higher gains make the jump to the left of point C, both of which are close to the desired yd=-3.5. This was verified and in each case an almost exact solution was found.
The above case studies and many more carried out and not reported here demonstrate the effectiveness of the proposed method in dealing with ill-condition cases, including severe nonlinearities and discontinuities.Figure 4
Case3: the reference trajectory yr(t), (a), actual trajectory y(t), (b), error e(t) (c), and inverse solution x(t) (d).
(a)
(b)
(c)
(d)


## 5.2. Example 2

We consider the complex two-input, two-output system described by(37)y1=sin⁡xcos⁡xy2=x12+2x22-x12x2;-3<x1,x2<3.The plots of the above functions are shown in Figure 5. For a desired set of (y1d,y2d), there can be multiple inverse solutions, and arriving at a particular solution depends on the initial values.Figure 5
Three-dimensional views of the functions (a)y1=sin⁡xcos⁡x and (b) y2=x12+2x22-x12x2.
(a)
(b)Case 1.
Lety1d,y2d=(0.174,1.22), x10,x20=(0.2,-0.5), y10,y20=(-0.4688,0.5600), and τ=20,k=0.5, and use the truncated singular value method to remove one singular value when it becomes below 0.01. The responses are shown in Figure 6, and the steady-state value of the output is obtained exactly equal to the desired values, and the corresponding input vector is x1d,x2d=(1.2,0.5).
We have computed the inverse desiredx1d and x2d for all values of y1d in the range of −20 to 60 and y2d from −1 to 1. The plots of inverses x1d and x2d against y1d and y2d are shown in Figure 7.Figure 6
Case1, two-input, two-output system: references trajectoriesy1r(t), y2r(t) (a), actual trajectories y1(t), y2(t) (b), errors e1(t), e2(t) (c), and inverse solution x1(t), x2(t) (d).
(a)
(b)
(c)
(d)Figure 7
Three-dimensional views of the inverse functions. (a)x1d as a function of y1d and y2d; (b) x2d as a function of y1d and y2d.
(a)
(b)Case 2.
The purpose of this case is to demonstrate the optimization discussed in Section4. Consider the two-input, single-output system of second equation in (37); that is, (38)y=x12+2x22-x12x2.Let the desired value be yd=1.22 and the initial values of x10,x20=(0.2,-0.5), as before. The corresponding initial output is y0=0.56. The optimization criterion is chosen as the minimization of ψ(x)=1/2(x12+x22) while achieving the desired output. Thus, gx=-∂ψx/∂x=-x1x2, and the parameters are chosen as τ=10,k=0.5, and β=1. The trajectories are shown in Figure 8 and indicate a smooth transition to the final values of x1d,x2d=(0.272,-0.7387) with the ψx=0.619. It is noted that with the same initial values but without optimization ψx=10.23. Thus, the optimization has minimized the change in x1,x2 while achieving the desired value.Figure 8
Case2, two-input one-output system with optimization: the reference trajectory yr (a), actual trajectory y(t) (b), error e(t) (c), and inverse solutions x1(t), x2t (d).
(a)
(b)
(c)
(d)


## 6. Conclusions

A novel approach to solving a class of nonlinear inverse problems is proposed. The solution is based on using a feedback system approach which ensures the stability of the solution using the classical control theory. Thus, the convergence of the solution is achieved which guarantees close to zero steady-state error, and therefore highly accurate inverse solutions are obtained. Two methods of regularization have been investigated within the feedback system: one is based on truncation and the other utilizes the Tikhonov method. Both of these methods have been shown to produce inverse solutions with a desired accuracy. It is shown that the proposed feedback scheme can also reduce the noise to a desired level. The proposed method has a number of features and advantages. A main advantage is dealing with highly nonlinear and discontinuous inverse problems. Furthermore, all results obtained for linear inverse problems such as dealing with ill-conditioned cases and regularization can also be achieved with the proposed formulation as a special case. The extension of the method to other classes of nonlinear inverse problems is currently being investigated.



---

*Source: 6843614-2017-05-28.xml*