Hi Saeb,

I may not be answering your question as you might expect, but just to clarify that in Mechanical Engineering, if you recall the Drag Force formula in Aerodynamics, then you can model the dynamics of a falling object.

During a free-fall, the net force acting on a falling object is the difference between the downward and upward force:

Fnet = Fdown − Fup

where Fdown is due to gravity and Fup is due to air resistance (drag). Using Newton's second law, F = ma:

m·v' = m·g − k·v2

where m is the mass of the falling object, v' = dv/dt, g is the acceleration due to gravity, and k is the drag proportionality constant determined by ½·ρ·A·CD. See the link for details.

https://en.wikipedia.org/wiki/Drag_equation

Simplifying the derivation from Newton's law and the ODE can be written in terms of differential rate of change (v'):

v' = g − (k/m)·v2

Your ODE looks like a kind of Bernoulli Differential Equation or Abel's Differential Equation of the first kind. Please check.

https://en.wikipedia.org/wiki/Abel_equation_of_the_first_kind

You can find Abel ODE representation of Friedmann equations (in physical cosmology) in this article:

Article Friedman versus Abel equations: A connection unraveled

Dear Yew-Chung. Hi, thank you for your concise and elaborate answer. Now, I am in possession of the answer to my query. Bests, Saeb.

I'm curious, what exactly are you interested in this 1st-order nonlinear ODE?

ẋ = x3 + u

@ Yew-Chung Chak

As part of my researchgate project, titled 'snap feedback control', I am investigating on global stabilization of nonlinear systems. The candidate snap feedback law; u=- k x exp(-sqrt(abs(x))), fails to stabilize the scalar control system; xdot=x^3+u, while for large initial conditions, the simulation becomes unstable. I want to know the physics background of the model; xdot=x^3+u, to get it revealed, what leads to the simulation instability?!. Through your answer, I concluded, the system xdot=x^3+u, is probably only a purely mathematical system, not a physical (as my own judgement). The candidate snap feedback law, globally stabilizes a wide range of typical nonlinear control systems.

Hi Saeb,

In the real-world situations, the operating range of the physical system is limited. Although your proposed controller is not a global stabilizer, if the operating range of {ẋ = x3 + u} is [−1, 1], then the nonlinear system can be stabilized when the control gain k in your "Snap Feedback Controller" is selected such that

k > e1 ≈ 2.71828.

This critical gain value is obtained by solving

x3 − k·x·exp(-√|x|) = 0, where x = 1.

You can try to simulate this on https://octave-online.net/

function xdot = f(x, t)

xdot = zeros(1, 1);

k = 3;

u = - k*x(1)*exp(-sqrt(abs(x(1))));

xdot(1) = (x(1))^3 + u;

endfunction

x0 = [1],

t = linspace(0, 10, 200);

y = lsode("f", x0, t);

plot(t, y)

@Yew-Chung. Thank you friend for the illustration.

k > e1 ≈ 2.71828: I think this condition holds for xdot

Hi Saeb,

The condition k > 2.71828 holds for −1 ≤ x ≤ 1. Change the initial condition x0 = −1 and run the simulation. See Figure 1.

To verify the condition, solve x^2 − exp(1)·exp(-√|x|) = 0. See Figure 2.

To view the magnitude of xdot, just plot x3 − exp(1)·x·exp(-√|x|). See Figure 3.

To find out the extremum points of xdot, solve d/dx {x^3 − exp(1)·x·exp(-√|x|)} = 0. See Figures 4 and 5.

@Yew-Chung. Thank you for illustration of the problem. It is now totally comprehensible.

Hi Saeb,

I just want to clarify one more thing. If you want to design a globally stabilizing 'Snap Feedback Controller' {u = −k·x·exp(-√|x|)}, then the gain k should NOT be a constant and must be expressed as a function of x:

k = f(x).

There are many candidate functions for k. Here are just some:

k1 = √2·exp(x2)

k2 = (2·x2)(2·|x| + 1)

k3 = (2·x2)/exp(-√|x|)

The comparison between the controllers of different candidate functions is plotted on Fig. 1.

Note that the first one (k1) will result in an extremely large control output for a relatively large initial value of x0. Theoretically, it works but it is not suitable for practical implementation.

The second candidate (k2) is reasonably acceptable, but slowly converges to zero upon reaching ±0.5. However, the issue of slow response can be easily fixed by adding a proportional term: u2 = − k2·x·exp(-√|x|) − x. See Fig. 2 for the time response without the P-term and Fig. 3 with the P-term.

The final candidate (k3) is for observation purposes because u3 = − k3·x·exp(-√|x|) = −2·x3.

@Yew-Chung Chak. Marvelous analysis. Thanks.

@Yew-Chung Chak. Please let me know, why you said, k3 is suitable for observation. I repeat again, why for 'observation' purpose, and probably not for stabilization.

Because k3 leads to u = –2*x³, which doesn't show any feature of your Snap Feedback Controller and it has the same structure as the nonlinearity term, x³.

In other words, applying k3 has the same effect of using the common nonlinear stabilizing controller u = –2*x³, which is simpler. The meaning of observation in this context is for comparison purposes.

@Yew-Chung Chak . That's it. I wrongly thought you meant k3 is for observer design.

Yup, forget about it. I should have picked a more appropriate word to describe that as a comparison between the common nonlinear control law and the other candidate functions in your Snap Feedback Controller.

In the previous case, it has been shown that your Snap Feedback Controller can handle unbounded polynomial function such as x3. Now, you can work out some candidate functions k = f(x) to deal with other nonlinear systems that contain unbounded transcendental function such as

ẋ = exp(x) + u

ẋ = exp(x2) + u

because your Snap Feedback Controller is inherently exponential. Moreover, it is not easy to control such systems using polynomial-based control laws.

Of course, for local stability, you get a truncated Taylor series expansion of the transcendental function and then stabilize each terms using polynomial-based control law, but the point is to showcase your Snap Feedback Controller.

Following

Dear Yew-Chung Chak

The siso control system; xdot=-x^3+u, is called "the hyper-sensitive system", and is an important benchmark for testing optimal control methodologies, particularly numerical methods. But, xdot=x^3+u, has another story. Kind regards, Saeb.

Hi Saeb AmirAhmadi Chomachar,

Have you published your Snap Feedback Controller?

Hi, Yew-Chung Chak

Nope, I've not published on that topic. I am trying to publish on another research topic, related to robust global output feedback stabilization of nonlinear systems. The preliminary research is available at my project: https://www.researchgate.net/project/hyper-exponential-control-theory-HECT

Regards, Saeb.

Hi Saeb AmirAhmadi Chomachar,

Could you demonstrate how you would apply HECT to the control problem of the inverted pendulum on a cart?

To simplify the problem, a linearized model (based on small angle assumption) and the following parameters are given:

M = 1.0 kg

m = 0.1 kg

l = 1.0 m

g = 10 m/s2

Dear Yew-Chung Chak, Thank you for the interest and query. I've provided my solution as an image attached to this answer, kindly for your notice. Regards, Saeb.

Dear Yew-Chung Chak If there is no need to control x, but only theta is to be controlled, please notice my revamped strategy attached as an image to this answer. Kind Regards, Saeb.

If you control x simultaneously with theta (the first strategy), then instability and divergence of the control process is plausible, probably due to specific dynamics of the inverted pendulum, as I experienced it myself for nonlinear case simulation through an optimal strategy. Moreover, angle control (the second strategy) seems stable and good, through an optimal control strategy. I should add it that, I've not simulated HECT for your benchmark linear inverted pendulum, based on the solutions I've provided as images attached to this query. If you made the computer simulation, please let me know your viewpoints. Regards, Saeb.

Hi Saeb AmirAhmadi Chomachar,

Thanks for your idea. I don't think I can solve the equation algebraically without incurring the singularity in the control law. But I think it is possible to solve it more conveniently in this form (see image above).

Dear Yew-Chung Chak , thank you. Your idea seems great. I didn't know that the control-law will be singular. Anyway, I think your adjustment is reasonably justified and sounds interesting. Regards, Saeb.

Hi Saeb AmirAhmadi Chomachar,

Let's take a simple double integrator, θ'' = u as example. I mean, solving the equation algebraically 'by force' will lead to an unwanted situation known as 'division by zero' when the angular rate θ' = 0.

Dear Yew-Chung Chak

I think this is due to the glitch of the Newtonian Dynamics, which assumes F=ma, and M=Jθ'', or for m=1, J=1, it assumes θ'' = u. So the Newtonian model might not give a physically controllable system. I've doing a project on critically modified Newtonian dynamics (CMOND) at:

https://www.researchgate.net/project/critically-modified-Newtonian-dynamics-cmond

where I assume F=ma+t*mdot*a+t*m*adot. Then your double integrator model becomes θ''+tθ'''=u, which can dodges the singularity when the angular rate θ' = 0. I am not sure about that, but it's my own idea about your query.

Hi Saeb AmirAhmadi Chomachar,

I'm not a dynamicist, so I can't tell. But so far, the Newtonian Dynamics has served me well in my ordinary modeling of motion systems. With CMOND, do you imply that the singularity can be avoided by solving this modified equation? Rather than putting the focus solely on the left-hand-side, perhaps there is an underlying theoretical framework of your HECT on the right-hand-side, which is not fully uncovered.

Dear Yew-Chung Chak Your analysis capacity, is admirable. That's amazing. I've not thought about the challenge of θ' = 0. So it's a bit troublesome.But still I think, your sliding surface strategy works well for the system θ'' = u, even if θ' = 0, and θ=0. The challenge of stabilization of an inverted pendulum or a swinging pendulum is presented in the paper attached as a pdf.

There you would find that, although the energy of the pendulum could be kept constant or slowly varying, however, it seems θ is not stabilizable for nonlinear pendulum. Also for nonholonomic systems, as the Brockett integrator, the x3 coordinate usually is hard to get stabilized, so x3, should be damped faster than x1 and x2, to avoid such a challenge, otherwise x3 would damp very slowly and almost get constant profile. Moreover, Brockett proved that the system is not continuously feedback stabilizable, and probably the continuous feedback-law would be always singular at origin. I think the same hurdle, persists for inverted pendulum. However, I think your sliding surface hypothesis is still reliable and without the singularity challenge. I am interested to know your viewpoints. Thank you, Saeb.

A simple idea is to assume hyper-exponential damping of θ as:

θ' =- θ e^(e^-t)).

Differentiate once to get:

θ'' = u = - θ' e^(e^-t))-θ d/dt(e^(e^-t))= θ e^(e^-t)* e^(e^-t) - θ d/dt(e^(e^-t)).

Hence, u=U(θ,t), is easily found. This is a time-varying feedback, only in θ and t. That's it.

Hi Saeb AmirAhmadi Chomachar,

I think the swinging pendulum model in Chung (1995) can be used to describe the simple motion of an overhead crane with a rigid hoist cable. Due to gravity, the payload swing angle of the pendulum is self-stabilizing that will always end up in the equilibrium state (θ = 0°) no matter what the initial states are. Thus, only the cart position x needs to be controlled and it is rather easy.

In the practical situation, the motion of the overhead crane may induce the payload oscillation (swaying). Depending on the payload mass and the trolley speed, the payload and hoist may exhibit double-pendulum behaviors.

https://thewikihow.com/video_5FfrKoToA0c

On the other hand, the inverted pendulum on a cart can be compared to the commercialized Segway two-wheeled, self-balancing personal transporter. The system does not have singularity issues.

Dear Yew-Chung Chak Thank you for the elaboration. I've got new insight on the problem in terms of its real practice. Notice that, I've added an image attached to this answer, regarding application of HECT to double-integrator system which you already posed as a question. I am looking forward to have your viewpoints.

Dear Yew-Chung Chak Although the swinging pendulum (in overhead-crane example) is self stabilizing, but to move the cart and pendulum to a new position, could might have been dramatically affected by the coupled dynamics of pendulum, payload and base support. I think it is not very easy to control the motion in nonlinear mode. As you might have noticed, in Chung (1995), it has been mentioned below the title that:

An experimental investigation into the dynamics and control of a swinging pendulum on a cart leads to new theoretical insights into the control of nonlinear systems about a periodic orbit.

I think this is an important assertion, and there are some challenges left in research for that question.

Hi Saeb AmirAhmadi Chomachar,

Thanks for showing your steps in getting the control law for the double integrator. Firstly, you can see the red-circled states (θ), those are 'stiffness' parameters. A 'damping' parameter (such as θ') is needed to 'dissipate' the kinetic energy of the double integrator system, and it does not appear in your control law. Thus, θ won't converge to zero.

Secondly, those blue-circled things (η2 and η) as you defined them, are positive coefficients. Hence, the 'positivity' of your proposed control law, u(θ, t) implies that energies are being 'accumulated' in the double integrator, leading to the divergence of θ, until the hyper-exponential term e^(e^-t) dies out.

Dear Yew-Chung Chak Thank you for the comments. You may use the law (2) as the control-law which includes the damping term θ' . I've replaced θ' from (1) into (2), to get (3), but this is not necessary. Hence, view the right-hand-side of (2) as the feedback-law u; which contains θ, θ' , t, and do not expand θ' in terms of θ.

Hi Saeb AmirAhmadi Chomachar,

Perhaps you overlooked the second paragraph of my previous post. Another issue is the sign of the stiffness term. Positive sign implies the accumulation of energy. With the dissipation of kinetic energy and accumulation of potential energy simultaneously, θ will oscillate. Hence, if the sign of the stiffness term in Eq. (2) is changed to negative, the system will be stable.

Dear Yew-Chung Chak

I had read your whole comment. Notice that, the second term in right-hand-side of (2) is equal to: - θ' exp(-t), then the control-law would be:

u= - θ' (η*exp(exp(-t))+exp(-t)).

or in other words:

u= - θ' (-θ'/θ+exp(-t)).

In this way both θ and θ' which are the stiffness and the damping terms, would appear in the control-law; u. But I think for θ=0, the control-law becomes singular. The challenge is getting a bit weird and interesting. It does worth to be discussed more. I am looking forward to have your viewpoints. Please if you find free time also check my manipulation for u, presented in this recent answer, in the lines above. Thank you.

Hi Saeb AmirAhmadi Chomachar,

Once you have proposed control law (for example, Eq. (2), no back-substitution), I think the appropriate way is to perform a stability analysis. Even though it is a 2nd-order nonlinear ODE, try running a stability test to check if it complies with the Routh–Hurwitz stability criterion.

The closed-loop double integrator system will have the following characteristic polynomial,

P(s) = s2 + a1·s + a0

where

a1 = η·e^(e^(−t))

a0 = − η·e^(e^(−t) − t)

For simplicity of exposition, perform the stability analysis at t = 0.

Dear Yew-Chung Chak

My latest control-law for the double-integrator system is:

u= - θ' (-θ'/θ+Ln(-θ'/θ)).

It is of pure feedback type. I will discuss it more with you later.

Hi Saeb AmirAhmadi Chomachar,

The solution for this 2nd-order nonlinear ODE, θ'' = − θ'·(− θ'/θ + log(− θ'/θ)), is given by θ(t) = c2*e^{Ei(e^{c1 − t})}, where Ei(·) is the exponential integral function. Although the solution converges, the real solution (without complex numbers) is only valid for (θ > 0; θ' < 0), and (θ < 0; θ' > 0). Continue to develop the theoretical framework of your Hyper-Exponential energy dissipation theory.

Dear Yew-Chung Chak

As it seems, the challenge is going to be fruitful. I've left you a direct message to your Research Gate inbox. Regards, Saeb.

Dear Yew-Chung Chak

Notice that, in sliding mode control (smc), there is an important concept called "reaching-law". As I know, the sliding line or surface is in the second and fourth quadrants of the state-space coordinate system. The reaching-law allows the phase-space trajectory to come from first or third quadrants, to second and fourth quadrants, and then to attract it to the sliding-line or sliding-surface, by a proper sliding-mode control-law. Hence, the challenge; θ' θ

u =− θ'·(− θ'/θ + log(− θ'/ η θ))

Hi Saeb AmirAhmadi Chomachar ,

I think you don't need to 'invoke' the sliding mode action yet. You can manipulate the 'minus sign' (−) so that your proposed control law is still valid when (θ > 0; θ' > 0), and (θ < 0; θ' < 0). Just replace all (−) signs with the sign function, sign(q), where q = θ'/θ, and you'll get:

u = sign(q)·θ'·{sign(q)·q + log[sign(q)·q]},

which can be simplified to

u = sign(q)·θ'·{|q| + log(|q|)}.

Although the modified control law looks too good to be true – there is a catch! It only works so long as θ'/θ < 1. This condition implies that approximately 75% of all operating points in the phase space can display asymptotic convergence (reaching a small neighborhood of zero). Look at this test simulation where θ(0) = π/20 ≈ 0.15708 rad, and θ'(0) = 0.157 rad/s.

Hope you gain closure with this new finding.

Dear Yew-Chung Chak

Thank you. Marvellous analysis. That was a nice job. Please also hint me why

θ'/θ < 1, is a pre-condition for the control-law.

Hi Saeb AmirAhmadi Chomachar,

The precondition θ'/θ < 1 is implied in this Vector Field (see image). If it is violated, θ will diverge (as you can see the quivers pointing away from the origin). Technically, you should perform a rigorous stability analysis to show this condition. By the way, do you think the origin is a stable node or a saddle point?

Dear Yew-Chung Chak

I think, whether the origin is a stable node or an unstable node, depends on initial conditions of θ, θ'. For some initial conditions it is a stable node, and for some other initial conditions it is an unstable node. It seems in its totality, origin behaves like a saddle point, stable for some initial conditions, and unstable for some other initial conditions, as you already said through this thread discussion:

"This condition implies that approximately 75% of all operating points in the phase space can display asymptotic convergence".

I think when you adjusted the control-law by assuming -1=sign(q) =sign(θ'/θ), then the phase-plane portrait (the vector-field image you illustrated) will change and the trajectories will tend to the origin in the first and third quadrant (against their current behavior as they leave away from the origin). It seems origin behaves like a saddle point mixed with limit cycles.

Meanwhile, I still do not claim that my hypothesis is correct. Looking forward, to know your comments on saddle-point behavior or stable node.

Hi Saeb AmirAhmadi Chomachar,

I simply don't know. The graphic on the vector field resembles a saddle-point, but 75% of the operating points in the phase space converge to a small neighborhood of the origin. The saddle points are always unstable. Although the quivers in the 1st and 3rd quadrants are leaving the origin, 50% of those quivers return to the 2nd and 4th quadrants, converging to the origin. Initial simulations showed they are stable.

You may perform some tests to find out if the origin satisfies the saddle point criterion. Because your control law is one-of-a-kind and unconventional, further theoretical framework is needed. If imaginable, find any physically-interpretable dynamical system that can be described by such ODE.

Dear Yew-Chung Chak

Thank you. It requires some careful theoretical investigation, and at present, I've gotten no more idea about that, and I don't simply know the answer alike you.

Dear Yew-Chung Chak Hi,

I've gotten an interesting result by assuming a time-varying η as; η= η(t)=e^e^(-t)

or simply; η=exp(exp(-t)). Assuming η to be constant, was the glitch of the previous solution. The control-law is changed from pure feedback type, to time-varying feedback, however, the good thing is that, the whole state space points will be globally stabilized to the origin in finite time, for any initial condition. Please notice the image attached.

Hi Saeb AmirAhmadi Chomachar

It's good to have discovered something. I couldn't find the stability analysis in your image. Could you show the basis of the global stability? Assuming k = 1 and σ = 1, try analyzing the following theoretical initial conditions {θ(0) = 1, θ'(0) = 2.8}.

How did you come up with the idea of using the logarithmic function? I have wanted to ask you that for so long. Logarithmic functions can produce complex numbers (z = x + i·y) and tend to minus infinity (−∞) when the argument approaches zero.

Dear Yew-Chung Chak

Thank you for your comment.

Please for the theoretical initial conditions {θ(0) = 1, θ'(0) = 2.8}, assume: k=10, and σ=1, to get:

η=10*exp(exp(-t)).

Matlab Mfile code:

function dy = rigidw(t,y)

dy = zeros(2,1);

dy(1)=y(2);

dy(2)=sign(y(2)/y(1))*y(2)*(sign(y(2)/y(1))*y(2)/y(1)+log(sign(y(2)/y(1))*y(2)/y(1)/(10*exp(exp(-t)))));

[t,y] = ode23s(@rigidw,[0 100],[1,2.8]);

Notice that ode45 seems with slower simulation run-time and almost unstable. I used ode23s.

The reason for such η=10*exp(exp(-t)), is in a very subtle point, concealed in the operation of stabilization. If you want to get an idea, look at my previous comment and the image, where I have noted at the end of the image, the point that: η=-θ'(0)/(θ(0)e). Hence, if η is constant for example η=1, then for θ(0) = 1, the only allowed initial condition for theta-rate is: θ'(0) = -e. Hence, for constant η, if you set one of the θ(0) or θ'(0), the other initial condition should be a certain prescribed constant with NO freedom to vary, and this is a philosophical glitch of the stabilization process which deteriorates the feedback-stabilization procedure. It means θ(0) and θ'(0) are not independent, but interconnected, for a constant η.

Hence, if you allow η to change with time, as: η=10*exp(exp(-t)),

then since η has an interval of variation, then the condition: η=-θ'(0)/(θ(0)e),

would remain logical. For different initial conditions, you should vary η with time, to cover a larger part of the state space and initial conditions. To vary η dynamics with time, to cover such bothersome initial conditions you mentioned as:

{θ(0) = 1, θ'(0) = 2.8}, you should set a suitable σ and k, as I set for an example above; k=10 and σ=1.

If you want to know how I've reached: η=k*exp(exp(-σ *t)),

it is based on the hyper-exponential damping of θ in Eq. (1) of the older image above this thread, where I assumed: θ'=-η*θ*exp(exp(-t)),

and hence, the option: η=k*exp(exp(-σ *t)), does not ruin the hyper-exponential dissipation of θ in Eq. (1).

I am looking forward to have your viewpoints. Regards, Saeb.

Hi Saeb AmirAhmadi Chomachar

Thanks for your explanation. I think I get what you want to do with η in order to ensure that the argument of logarithmic function is less than 1. This approach requires a priori assumption about the initial conditions. In practice, the operating ranges for {θ, θ'} are generally known. Is it possible to design a robust η(t) that covers -1 ≤ θ ≤ 1 and -1 ≤ θ' ≤ 1?

By the way, I'm interested to know how you came up with the idea of logarithmic function (as circled on the image above) in the control law. The complex numbers (z = x + i·y) and minus infinity (−∞) issues restrict the logarithmic control law from functioning at certain regions of the phase space.

Second matter: how did you come up with the idea of quotient θ'/θ in the control law? Conventional wisdom tells the control engineers to avoid or prevent singularity. If the logarithmic control law commands at infinite force/torque, the physical actuator will behave problematically.

Dear Yew-Chung Chak

The term you circled with red color, is f=Ln(d), where f=e^(-t), and since:

θ' = -η θ e^e^(-t), was the primary assumption,

then; f=Ln(-θ'/(η θ))=e^(-t), where d=-θ'/(η θ). In this way, I removed the highly time-varying term, from the right-hand-side of Eq. (2), in the image posted already in older comments above this thread discussion. This changes the highly time-varying control-law to a semi-feedback control-law, I mean a time-varying feedback-law, with less dependency on time.

Please tell me why did you say:

If the logarithmic control law commands at infinite force/torque, the physical actuator will behave problematically. I think by this assertion, you meant when θ approaches infinity, then the term inside logarithm becomes zero, and logarithm becomes minus infinity ( −∞ ). Moreover, notice that, usually θ and θ' kick off at a known bounded initial condition and approach zero with the similar order of magnitude, at steady state, hence logarithm, would usually approach Log(ε/ ε =1) or zero, and it does not have the singularity issue in practice. However, I think robustness issue is still meddling when θ becomes very large, and θ' is finite, and hence the logarithmic term becomes minus infinity ( −∞ ). But still I think such a feed-back law is worthwhile, while other methodologies presented in the existing literature, to control double-integrator system, suffer from too many drawbacks, already known.

Please, also tell me what is the physical interpretation of the quotient θ'/θ , as you've already called it conventional damping and stiffness. Do you mean the terms θ, θ' are like damping and stiffness terms in conventional second order systems?

Although I first proposed k=20, σ=10, but notice that, even for k=10, σ=1, the stabilization process works well for the theoretical initial conditions {θ(0) = 1, θ'(0) = 2.8}, which you had already asked.

Hi Saeb AmirAhmadi Chomachar

Thanks for your explanation on the log term. C·θ' is the damping term and K·θ is the stiffness term. They are commonly used in the analogous mass-spring-damper 2nd-order systems by many mechanical engineers and materials scientists.

Regarding the quotient θ'/θ, if θ = 0, the quotient will result in division by zero. When θ rapidly approaches 0, the quotient tends to get large, as you can see the spike in the figure. Saturation will occur if the spike hits the output limit of the actuator.

The more important question is: Given the operating ranges for {θ, θ'}, how would you design a stabilizing η(t) that covers -1 ≤ θ ≤ 1 and -1 ≤ θ' ≤ 1?

Dear Yew-Chung Chak

Thanks for your query.

Since, the feedback-stabilization of double-integrator system, is sought for quite a long time in the control systems society, with little insight on the problem presented so far, I think the spike in control magnitude is not very discouraging, particularly through a theoretical scope of nonlinear control theory, where even relaxed control signals (highly-chattering), and also discontinuous control signals, are still theoretically worthwhile, and even practically they have some applications, although they seem dysfunctional through implementation scope.

Your question: Given the operating ranges for {θ, θ'}, how would you design a stabilizing η(t) that covers -1 ≤ θ ≤ 1 and -1 ≤ θ' ≤ 1?

It sounds a good question. I think although for the case; -1 ≤ θ ≤ 1 and -1 ≤ θ' ≤ 1, the situation is worsened particularly for θ' near zero and θ finite, or vice versa, which yield the logarithmic term as infinite, seen ultimately as a spike in control signal, but notice that as I know, Log(x) is Lipschitz for large x (x approaches ∞) and I think this fact ameliorates the challenge of unbounded control signal, at least theoretically. Even for small x, the term Log(x) is minus infinity or -∞, but I think in real-time dynamics, such a condition does rarely happen when we assume the dissipation rate and order of magnitude of θ and θ' in real-time. If this is not the case, then I think, we should turn for critically modified Newtonian dynamics (cmond) in my project:

https://www.researchgate.net/project/critically-modified-Newtonian-dynamics-cmond

which accounts for jerk variation, to get new insight on second order dynamics model of the double-integrator system. I think it is due to the Newtonian model inaccuracy that we reach such singularity of the control-law. If the model is updated based on modified Newtonian dynamics (cmond), then I guess the challenge of singular control and as well initial conditions partial coverage of state-space, would be apparently bypassed but I should exercise it to get certain, I am not sure.

If you have other novel ideas, I am humbly looking forward to know about them. Your comment is kindly welcome. Best, Saeb.

Hi Saeb AmirAhmadi Chomachar

I'm just wondering how you would design one generalized η(t) that stabilizes the entire phase space within the region bounded by -1 ≤ θ ≤ 1 and -1 ≤ θ' ≤ 1. The task can be onerous, if the operator has to re-tune η(t) for each new initial conditions, that violate θ'/θ < 1, in the operating system.

In this practical situation where θ' and θ are finite, what will be the k in η(t) = k·exp(exp(−σ·t))? Perhaps the answer can be satisfactorily found when you fully develop the critically modified Newtonian dynamics approach.

Dear Yew-Chung Chak

Thank you for the query.

Please explain more, what is the problem with θ'/θ < 1 ?

I simulated the control process with some initial conditions in the region θ'/θ < 1, while I set; k=1 and σ =1, and there was no problem with simulation. Notice that when k=1, then 1< η(t) -1/e, which hyper-conditionally satisfies the assumption:

1< η(t)

Hi Saeb AmirAhmadi Chomachar

The instability condition θ'/θ ≥ 1 is observed in the vector field (setting η = 1) as shown and explained previously. If θ'(0)/θ(0) ≥ 1, I know you can fix it by introducing η(t) so that θ'/(η·θ) < 1.

Both initial conditions {θ(0) = 1; θ'(0) = 0.5} and {θ(0) = 0.8; θ'(0) = −0.5} satisfy θ'/θ < 1, hence η(t) = 1*exp(exp(−t)) should work well. However, does the same η(t) work for the new initial condition {θ(0) = 0.36; θ'(0) = 1} as θ'(0)/θ(0) > 1 and bounded by −1 ≤ θ ≤ 1 and −1 ≤ θ' ≤ 1? Try running ode23s.

Dear Yew-Chung Chak

Thank you for the comment.

Beware that, for time-varying systems (with time-varying terms in the right-hand-side of the ode), as η(t) in the case of applying HECT to double-integrator system, the vector field plot of (θ, θ') cannot be reasonably interpreted, in my own viewpoint and insight.

I mistakenly thought you meant the challenge of initial conditions resting in the region θ'/θ < 1, whereas you meant the initial conditions which violate θ'/θ

Hi Saeb AmirAhmadi Chomachar

I think you have unknowingly discovered another 'gem'. I didn't realize that when η > 1, the instability condition is also changed. Hence, θ'/(η·θ) ≥ 1 no longer applies, as shown in the vector field. You have changed the dynamical properties into a 'magnetic field'-like attraction to the origin (see image). Bravo!

Dear Yew-Chung Chak,

Thank you. Would you please explain more about the 'gem' which I have unknowingly discovered ?

and also about your comment: "You have changed the dynamical properties into a 'magnetic field'-like attraction to the origin".

Please explain more about the new discovery. What is the important implication of the then discovery (not known by me) for the control process if η =2.7183.

Is the challenge of the initial conditions θ'/θ ≥ 1, not covered in the control process, bypassed and resolved in the way, and how? What would be the revamped η=η(t), and what about the values for k and σ. Please lay out more information about the discovery. I do not know, what you are talking about (?).

I am earnestly looking to have your reply. Regards, Saeb.

Hi Saeb AmirAhmadi Chomachar

It has been thought (by me only) that the condition θ'/(η·θ) < 1 must be satisfied (according the vector field). However, you introduced η(t) = 1*exp(exp(−t)) or η(0) ≈ 2.7183, which tears down the boundary θ'/θ = 1, and wraps the quivers back to the origin (like the magnetic field). See the previous image. Although I couldn't understand why it works, you can test {θ(0) = 0.0036; θ'(0) = 1} with a constant η = 2.7183.

Dear Yew-Chung Chak ,

Thank you. I simulated the Matlab Mfile:

function dy = rigidw(t,y)

dy = zeros(2,1);

dy(1)=y(2);

dy(2)=sign(y(2)/y(1))*y(2)*(sign(y(2)/y(1))*y(2)/y(1)+log(sign(y(2)/y(1))*y(2)/y(1)/(2.7183)));

with 23s ode solver:

[t,y] = ode23s(@rigidw,[0 1],[0.0036,1]);

and also with ode 45 solver;

[t,y] = ode45(@rigidw,[0 1],[0.0036,1]);

for the initial conditions you suggested as: θ(0) = 0.0036; θ'(0) = 1 ,

and for constant η = 2.7183.

But unfortunately both ode solvers gave unstable simulation result.

Would you please check more accurately the initial conditions and η = 2.7183, if there is something missed in the data you provided. I also suspect a buggy Matlab simulation, but I need your help with that to get sure about the error.

Moreover, as we already discussed η=θ'(0)/(θ(0)e), so how do you want to satisfy such condition, with your proposed data which gives η=θ'(0)/(θ(0)e)=10000/36/e , which contradicts your assumption that η = 2.7183 (?!).

Regards, Saeb.

Hi Saeb AmirAhmadi Chomachar

The quotient should be consistent.

η = 1·exp(exp(−0)); % η ≈ 2.7183

q = θ'/(η·θ);

u = sign(q)·θ'·(abs(q) + log(abs(q)));

It is unnecessary to make η = θ'(0)/(θ(0)·e) because of the 'gem (η > 1)' you unknowingly discovered. This condition θ'/(η·θ) < 1 no longer applies. Thought you already see the vector field.

Dear Yew-Chung Chak,

Thank you. I got the good simulation result. Please share with me more ideas about this accomplishment. I am looking forward to have your supplementary comments and viewpoints. I am obliged to you, Saeb.

Dear Yew-Chung Chak

Hi, thank you for your help. I am trying to get the quiver plot you illustrated as an image attached to your answer. I've used the Matlab syntax:

[x,y] = meshgrid (-1.5:0.2:1.5);

Z=sign(y./x/2.7183).*y.*(abs(y./x/2.7183)+log(abs(y./x/2.7183)))

[u,v] = gradient(Z);

h = quiver (x,y,u,v);

title ("quiver plot");

Where I assumed x= θ and y= θ', and Z= θ"=f(x,y).

However, I think your plot is more pellucid and the arrows size in your quiver plot are uniform. Would you please help me, what is the issue with my Matlab syntax which renders a different plot from that of you? Did you write your code in Matlab software environment? Is the assumption [u,v] = gradient(Z); in my syntax correct? Is your quiver plot a gradient field? Regards, Saeb.

Hi Saeb AmirAhmadi Chomachar

I ran the simulation on Octave Online. I just specified the velocity vectors u and v for x' and y', respectively. The size of the quivers are re-scaled for presentation.

[x, y] = meshgrid(-1.5:0.15:1.5, -1.5:0.15:1.5);

q = y./(1*exp(exp(-0))*x);

u = y;

v = sign(q).*y.*(abs(q) + log(abs(q)));

quiver(x, y, u, v, 'AutoScaleFactor', 0.1, 'MaxHeadSize', 0.08)

axis equal

axis([-1.5 1.5 -1.5 1.5])

xlabel('\theta')

ylabel('\theta''')

title('Vector field with \eta = 2.7183')

Dear Yew-Chung Chak,

Using Matlab, I got the plot through a larger window. Notice the image attached.

Thank you. I am obliged to you. Saeb.

My syntax duplicate from your Octave Online code:

[x,y] = meshgrid (-5:0.4:5);

u=y;

v=sign(y./x/2.7183).*y.*(abs(y./x/2.7183)+log(abs(y./x/2.7183)));

h = quiver(x, y, u, v, 'AutoScaleFactor', 10, 'MaxHeadSize', 20);

title ("quiver plot");

xlabel('\theta')

ylabel('\theta''')

global feedback stabilization of double integrator system

Dear Yew-Chung Chak

Hi,

Please let me know your idea about the quiver-plots for the optimal feedback stabilization of double-integrator system, attached to this comment. See the project:

https://www.researchgate.net/project/Analytic-Robust-Nonlinear-Optimal-Feedback-Control

Hi Saeb AmirAhmadi Chomachar

It looks a bit similar to the sliding mode. Have you tested your proposed controller for x(0) = 1 and x'(0) = 0?

Dear Yew-Chung Chak

Thank you for the comment.

As it is seen in the control-law, x2=0, is a point of singularity. The singularity at equilibrium point is so much common in control systems, that an equilibrium point is also called a singular-point in the terminology. But I should discuss it more with you.

For zero initial velocity, you are right. It gives indefinite simulation. I think a hybrid control-law could possibly resolve the challenge. At initial time; t0=0, move the object to a very small vicinity of its initial position, to make velocity a very small non-zero value. We can do this by applying a very small-gain linear feedback as:

uLinear=-k*sign(x1), where k could be small, and afterwards we apply the optimal feedback-law for the whole simulation time remaining. It means, the control-law could be adjusted to:

u=uLinear*(H(t)-H(t-c))+uOptimal*(H(t-c)), where H is the Heaviside step function and, c is a very small value corresponding to the time-interval when the Linear feedback-law is applied to slightly move the object, and make its velocity non-zero. For example c=0.1, c=.01, c= .03, or any other suitable value.

I tried to simulate the hybrid control-law for the double-integrator system, in Matlab Simulink software environment, but there was error notification due to some software-oriented reasons, NOT probably due to the control-law itself. I will try more, and let you know if I get the simulation bypass the error.

Dear Yew-Chung Chak

In the application of the HECT to double integrator, as I noticed, there was an η missed behind the first term in the parenthesis which is abs ( θ' /2.7183/ θ). I have attached an image for your notice. When I add the coefficient η=2.7183, before the abs value in the parenthesis, the time-domain simulation becomes unstable. Please comment on that, and rearrange the control-law with the missed coefficient η. I will be thankful.

Dear Yew-Chung Chak

regarding singularity of simulation for x'(0) = 0, at initial instance of control, we can apply the simplified control-law (u=-2*x1); for the case x2 approaches zero, and then apply the full feedback control-law. It means, at first, we initialize the system by a very minute movement, and then employ the fully nonlinear optimal feedback-law to bring the object to the origin.

Hi Saeb AmirAhmadi Chomachar

I'm unsure, but I think that you have breached your HECT stability condition by adding the 'missed' η = 2.7183. I don't know what the exact HECT stability condition is, but you have find it in order to determine the 'missed' stabilizing parameter. Only when it is stable, then you can rearrange the control law in a manner you deem appropriate.

Regarding the singularity issue in your Leibniz-based optimal nonlinear control law, are you suggesting that a piecewise controller, where the controller is defined by two sub-controllers and Each sub-controller is triggered depending on state x₂? If you do this, is the system considered truly optimal or sub-optimal?

Dear Yew-Chung Chak

It might not be truly optimal, but I think it is still a globally feedback stabilizer piecewise controller, and this is theoretically as well practically worthwhile.

Dear Yew-Chung Chak

Regarding missed η, in control-law, I should add it that, if we add the missed η to the control-law, the quiver-plot still as you interpreted is similar to magnetic-field like attraction to the origin, with a smaller slope. But, the time-domain simulation becomes unstable. What can we conclude from this observation? Thank you.

Hi Saeb AmirAhmadi Chomachar

I don't know what your initial conditions are, but the quiver plot is not a tool for demonstrating absolute stability. Your initial conditions maybe lying on the hidden instability regime, and the quivers are surrounding the hidden instability regime.

The quivers do not show every single points in the plot. Therefore, I still think that you might have breached your HECT stability condition with the added 'missed' η = 2.7183. Determine the stability condition is next important task to do. It is essential to note that a system that is stable for gain K = 1 may become unstable for a different gain K = 2.7183.

Dear Yew-Chung Chak

Thank you, for the comment. I should find some ideas to resolve the challenge. I would let you know if there were any breakthrough.

Dear Yew-Chung Chak

Regarding the Leibniz-based optimal control-law, I think practically we can assure initial zero-velocity never happens, although initial velocity could be very low as 10^(-100). Hence, apart from the theoretical challenge, the proposed optimal control-law is resourceful, if we set the initial velocity to a very small value. Please, see my recent query at the URL:

https://www.researchgate.net/post/Is_there_any_inert_object_in_our_world

Hi Saeb AmirAhmadi Chomachar

You are right that there might be no absolute zero velocity, but I think the state x₂ in your Leibniz-based optimal nonlinear control law, refers to the relative velocity. Hence, if an object appears stationary with respect to a reference frame, its state x₂ = 0 m/s or rad/s.

Dear Yew-Chung Chak

With regard to the optimal control law, please let me know,

how you interpret the situation illustrated in the image attached here, for a generic control practice. Is x2(0) simultaneously zero and non-zero at the initial time t0?.

Hi Saeb AmirAhmadi Chomachar

Is this some kind of a 'Paradox'? As far as I know, contradictory initial conditions cannot both be true in the same sense at the same time. Could it be a misprint in the illustration? Otherwise, I believe that the safest way to interpret the information is shown in the image above.

Dear Yew-Chung Chak

I think this is a paradox, as you illustrated by the hypothetical red line, connected to the x2 axis, as it becomes zero at a non-zero point. It means graphically, the vector at initial time says: dx1(t=0)=0, hence, dx1/dt(t=0)=x20=0, but through the phase-plane coordinate; x20 is not zero, and you showed this paradox by the red line connecting to the x2 axis. I am still not in full possession of answer to this question, however, I know that, in optimal control, the philosophy of initial conditions differs from conventional classic control. For concreteness, in optimal control-law, usually x(0) and x(T), denoting state values at the beginning and the end of operation, is explicitly incorporated in the optimal control-law; for example as u=x0*x^3-x(T), this is only an example, not necessarily meaningful. Whereas, in classical control, usually such initial conditions are not seen in the control-law. Also notice that, optimal control, particularly when numerical, is very sensitive to initial and end conditions (terminal set). The mathematical philosophy of initial and terminal conditions in optimal control is much deeper than to judge it at the first glance, and as I said, I still do not fully understand it.

Dear Yew-Chung Chak

I was just tempted to try your sliding-surface adjustment to HECT, and apply it to double integrator system. So I did it, and simulated in Matlab, while the illustrations are presented as attachments to this comment. It seems the system is globally hyper-exponentially stable (with any initial condition), but the control signal becomes non-zero in steady state.

I guess you might have observed this phenomenon yourself, before you had asked me to apply my HECT theory to double integrator system, at the midst of this thread discussion. Am I right? Have you already tried it? If yes, please let me know your comments about this phenomenon.

I had a glance on the 2018 paper:

"On hyper-exponential output-feedback stabilization of a double integrator by using artificial delay"

by: Efimov D, et.al.

and it seems the problem of hyper-exponential stabilization of double integrator is not so simple.

I've gotten some new ideas to apply HECT for double integrator system, based on cost minimization. I would let you know about that, if I get fruitful result.

Hi Saeb AmirAhmadi Chomachar

I didn't simulate this, but simple analysis

• u(t, x) = e^e^(−t)·{[e^(-t) − e^e^(−t)]·θ − 2·θ'}
• u(t, x) = f(t)·{h(t)·θ − 2·θ'}
• f(t) → 1 as t → ∞
• h(t) → −1 as t → ∞
• u(t, x) → 1·(− θ − 2·θ') as t → ∞

showed that the double integrator system eventually becomes

θ'' + 2·θ' + θ = 0 as t → ∞. I guess you know where θ'' + 2·θ' + θ = 0 will lead to.

Dear Yew-Chung Chak

I found the solution to the ode:

θ'' + 2·θ' + θ = 0 ,

as:

θ=c1.exp(-t)+c2.t.exp(-t),

and hence:

θ''=u = c1.exp(-t)-3.c2*exp(-t)+c2.t.exp(-t).

and I guess due to explicit time variable t inside u and as well θ, the control implementation of such a system is troublesome in practice.

But with regard to your comment:

"I guess you know where θ'' + 2·θ' + θ = 0 will lead to".

Is there any other issue this ode will lead to? So, please let me know.

Dear Yew-Chung Chak

I think, I've found the glitch. Please, notice the image attached.

Moreover, I have gotten some new ideas on simple application of HECT to double-integrator system, and as I'd already promised, I will let you know, if there were any fruitful result.

Hi Saeb AmirAhmadi Chomachar

{θ'' + 2·θ' + θ = 0} is a stable ODE that is converging towards (0, 0). Your claim that the control action is non-zero at steady state appears to contradict the finding.

θ' = − η·θ such that η > 0, is an exponentially stable 1st-order ODE. I don't know your purpose of the differentiation to obtain the unstable 2nd-order ODE θ'' = η2·θ. If you want to ensure that a 1st-order ODE θ' to have a stable 2nd-order ODE θ'', then use this one:

θ' = θ/(1 + t) − η·θ

Dear Yew-Chung Chak

I wanted to showcase it that, with direct differentiation of hyper-exponential equation, as we did it before above this thread discussion, to get θ'' =u, in terms of θ' , θ , t , we cannot guarantee the stability of the second-order system based on the hyper-exponential stability of the first-order system (ode). Hence, our derivation is not dynamically justified (admissible) although mathematically feasible. I only made the example θ' = − η·θ, and differentiated it again to show how we end in an unstable system. This is an example for the case θ' = − η·θ e^e^(-t), and it seems its second differentiation ends in an unstable double-integrator control system.

Another question which is not a red herring, why don't we select u=-2θ'-θ, to simply control the double-integrator system, as θ'' =u?. Is the main challenge for convergence time? Or linear feedback is not suitable and robust? What is the challenge of linear feedback to simply control the double-integrator system?

Hi Saeb AmirAhmadi Chomachar

I think we have to flex our brain muscles every now and then. You have proposed the unconventional HECT approach for the general 2nd-order systems, and you are thinking some ways to handle the singularity issue. Perhaps, modifying your HECT design approach would eliminate the division by θ' (velocity state). It is also important to note that from θ' = − η·θ to θ'' = η2·θ leads to instability.

The PD control, u = − kd·θ' − kp·θ is free from the singularity issue. I have been using it consistently for some special types of 2nd and 3rd-order systems and obtained good results. If there is a steady-state error, a little Integral action can be added. Moreover, the desired convergence time can be adjusted by tuning the PD gains. I think you have probably heard that the PID control architecture is the most widely used controllers.

'Re-connect' your emotive intuition when you first came up with your HECT:

• Why does it need to be a hyper-exponential function?
• What specific control problem that your intend to solve with HECT?
• What classes of systems can HECT be satisfactorily applied to?

Dear Yew-Chung Chak

Thanks for your explanations.

Fortunately, I've gotten new results on application of HECT to double integrator system. I need your advice to further adjust the procedure to get more satisfying simulation outcome.

I am looking forward to know your viewpoints and comments.

Please see the pdf attachment along with the illustration.

The simulation is performed for initial conditions:

θ(0)=1, θ'(0)=-1.

η= 2.7183

I've observed that for θ(0)=1, θ'(0)=1, if we adjust u with sign convention similar to what you already proposed above this thread discussion, as:-1=sign(η*θ* θ'),

then the simulation will be stable, but θ positively grows (and saturates) while θ' damps down.

Dear Yew-Chung Chak

For the case of initial conditions being both positive as:

θ(0)=0.3, θ'(0)=1,

the Matlab Simulink simulation with the control law adjusted with the sign convention as:

u=(-sign(η*x1*x2)*2*η*x2^2*exp(exp(-s*t))+s*sign(η*x1*x2)*2*η*x1*x2*exp(-s*t)*exp(exp(-s*t)))/(η*x1*exp(exp(-s*t))-x2);

with η= 2.7183, and s=1; be the exponential power coefficient,

where in u expression, x1=θ, x2=θ', and the plot is illustrated in the image attached, for your notice.

Hi Saeb AmirAhmadi Chomachar

I'm unsure, but for a double integrator, {θ'' = u}, I think there is a misconception in Eq. (4−H). Eq. (1−H) shows your desired hyper-exponential stability, described by the 1st-order derivative for θd. Eq. (3−H) is the derivative of Eq. (1−H).

Why would you equate θ" = θd" in Eq. (4−H) to obtain the control law (u) without verifying whether θd" is stable or not?

Dear Yew-Chung Chak

Thank you for your comment. My main question was with regard to the PDF attached to my query, which I guess you didn't notice it. Please read the PDF attachment and let me know your viewpoints. Kind Regards, Saeb.

Hi Saeb AmirAhmadi Chomachar

I thought you are open to feedback since they are part of the contents in the PDF. I was just trying to correct the misconception according to the control theory.

Regarding your four control laws, the first three control laws have singularity at θ = 0. The 4th control law has singularities at the locus that satisfy θ' − η·θ·e^e^(−t) = 0. I think all four control laws will lead to some saddle points. You need to verify that.

Dear Yew-Chung Chak

I missed to spot that the formulas were part of the PDF, I am sorry.

So, thank you for your comment. I think I should revise my idea of employing HECT for control of double-integrator system. As it seems, generalizing HECT to double-integrator system is a bit challenging, and not as simple as using HECT for system of first-order ode.

Hi Saeb AmirAhmadi Chomachar

Searching the correct methodology is very important. You picked a cost function J and you design a control law to minimize it. Hence, your control law will only guarantee the shortest distance in the trajectory from time 0 to t (see image for u1), but it does not guarantee the stabilization of the double integrator. In short, the missing ingredient in your HECT approach is the stability. So, your HECT has to be designed based on some stability theorems, not only the Leibniz integral rule.

Dear Yew-Chung Chak

Thank you for your thoughtful comments and answers regarding my queries. I will let you know if I find any breakthrough on the goal.

Moreover, please tell me what is the issue with the sliding surface plan and Lyapunov function as proposed a bit sooner above this thread discussion and attached here as an image.

You presented your analysis on stability of both the double-integrator system and the control signal, which is repeated here as:

Yew-Chung Chak wrote:

I didn't simulate this, but simple analysis

• u(t, x) = e^e^(−t)·{[e^(-t) − e^e^(−t)]·θ − 2·θ'}
• u(t, x) = f(t)·{h(t)·θ − 2·θ'}
• f(t) → 1 as t → ∞
• h(t) → −1 as t → ∞
• u(t, x) → 1·(− θ − 2·θ') as t → ∞

showed that the double integrator system eventually becomesθ'' + 2·θ' + θ = 0 as t → ∞. I guess you know where θ'' + 2·θ' + θ = 0 will lead to.

So is it OK to claim that, this control-law is a good alternative to control the double-integrator system, based on HECT?.

I guess you have no negative idea about this control-law. Am I right?

Thanks in advance.