A differential equation is an equation involving a function and its derivatives. Eg.
dx
-- = 4
dt
is a differential equation. Thus a differential equation is an equation whose solution is a function, or a family of functions. The above equation can be simply solved by integration, giving a value for x
where c is an arbitrary constant - ie the solution is a family of functions parameterised by c. An ordinary differential equation (ODE) is a differential equation in which differentiation is with respect to one variable only. A partial differential equation is a differential equation in which differentiation with respect to several variables may occur. Eg.
is a partial differential equation which describes many forms of wave motion. Robotics, for the most part, makes use of Ordinary Differential Equations, where differentiation takes place with respect to time.
Clearly, differential equations of the form
dx
-- = f(t)
dt
where f(t) is a known function of t can be solved by integration, which is not itself an operation that can always be performed analytically - sometimes integration gives rise to new functions. Thus solving differential equations is at least as hard as integration, and may sometimes have to be performed numerically. However for the control of robots (as in any other problem of Control Theory) we are not so much concerned to solve a differential equation as to choose a differential equation whose solution has certain desirable properties.
A differential equation is linear is linear if the unknown function(s) and its (their) derivatives occur linearly in the equation. Eg.
is linear, while
is not. Homogeneous linear differential equations have the property of superposition of solutions. Given an equation
??? an image missing here
an(t) \dxbydtn + a(t)n-1 \dxbydtnpred + ... +a0(t)x = 0
if x1(t) and x2(t) are solutions, and if c1 c2 are constants, then c1 x1(t) + c2 x2(t) is also a solution. (This is immediate from the fact that differentiation itself is a linear operator). The same property holds for simultaneous linear differential equations.
Linear differential equations with constant coefficients (ie where the an(t) are constants) are by far the most tractable kind of differential equations, and there is a rich body of knowledge about them. Therefore we try to cast our problems in a form in which these are applicable. This can be approached through suitable robot design (for example reducing friction and backlash). However, the differential equations for the dynamics of any robot with revolute joints will be non-linear, but for control purposes the technique of linearising the equations can be effective since the purpose of a controller is to ensure that the controlled system deviates from its intended behavior only by a small amount.
A simple example of linearising a differential equation can be seen in the case of a pendulum. The differential equation of motion is
whose exact solution requires elliptic functions. However the equation can be linearised, by approximating \sin(\theta) by \theta (measured in radians of course). In this case we get
m\dbydtsq{\theta} = \frac{m}{l}\theta$$
which is linear, with constant coefficients.
Consider a room, of heat capacity H, heated by a heater, of capacity W watts, We will assume that the temperature in the room has a uniform value T, and that heat is lost to the outside at a rate that depends on the temperature difference T-T0, of L watts per degree centigrade. If the heater is turned on then the differential equation for T is:-
while if it is turned off, then the ODE for T is:-
??? missing image $$ \dbydt{T} = \frac{ - L(T-T0)}{H} $$
Both of these are linear differential equations for T, with constant coefficients. The first equation can be written (with a positive denominator on the LHS):-
$$\frac{1}{a-T}\dbydt{T} = b$$
Integrating with respect to t we get:-
$$\int\frac{dT}{a-T} = \int b\, dt$$
So that:-
$$ -\ln(a-T) = bt+c$$
$$a - T = e^{-bt-c}$$
?????? picture missing here
$$T = a - de^{-bt}$$
where d=e^{-c}. That is T tends to a from below.
In the second case we have
$$\int\frac{dT}{T-a } = -\int b\, dt$$
with a positive denominator. Integrating, we get
$$ \ln(T-a) = -bt+c$$ $$ T = a + e^{-bt+c} = a + ke^{-bt} $$
So that T approaches a from above. Note that both the heating up and cooling down equations have the same form if we allow a positive or negative coefficient of e^{-bt}.
In the cooling case a=T0, so that the room temperature approaches T0 in the limit. In the heating case, let Tmax = T0+W/L, then T tends to Tmax in the limit.
If now we install a thermostat, which turns off the heater when the
T>T2 and turns it on when T
It is important to remember that a controller can only accomplish what is physically possible. A thermostat cannot by itself cool a room down below the outside temperature- you need an air conditioner.
We saw above that the solution of a what is a rather simple differential equation is rather hard work, and might despair of dealing with the complicated equations that can arise when we have multiple motors acting on a robot through flexible drive trains. In the case of linear ODE's with constant coefficients there is a very effective treatment through the following:-- Definition
If f is a function, then the Laplace transform
is defined by:-
$$F(s) = \laptrans{f(t)}$$
Clearly, there will be functions for which the Laplace transform is not defined because the integral does not converge, but such functions are not of significance from the point of view of robotics.
The reason Laplace transforms are so useful is in view of the following result, obtained by integration by parts:-
Now, provided e-stf(t) tends to zero as t increases unboundedly for some s this will evaluate to
This will indeed happen for any functions of engineering interest.
Thus, if we arrange, by adding suitable offsets if required, that all functions we apply the Laplace transform have the property that f(t)=0, we see that we can map the operation of differentiating a function into multiplying its Laplace transform by s.
Since, for practical purposes, functions are uniquely determined by their Laplace transform for the values of t>0, we can translate problems in differential equations into manipulations of Laplace transforms, which are often simple polynomials.
Laplace transforms can be determined by standard integration techniques. Eg.
{\cal L}(t) = \laptrans{t} = -\frac{1}{s}\intposx t \dfbydt{e^{-st}}\, dt$$
= -\frac{1}{s}\livb te^{-st} \rivbposx + \frac{1}{s}\intposx e^{-st}\,dt$$
= \frac{1}{s2}$$
In this section we tabulate the Laplace transforms of a number of functions. For control, important functions are the unit step:-
and the unit ramp:-
The following table can be obtained by standard integration techniques (the sin and cosine formulas require a double integration by parts).
Note that the Laplace transforms of all of these functions are rational functions of s. Moreover we can see that:-
(consider the differentiation rule and notice that the integral is zero if t=0).
So we see that for the whole vocabulary of elementary functions, and operations of differentiation and integration on them, the corresponding Laplace transforms are rational functions of s. Now any rational function over the reals can be expressed as a sum of partial fractions with linear or quadratic demoninators, or powers of the same. (Apply the fundamental theorem of algebra, and group complex conjugate pairs of roots together.) Thus we begin to see that it may be possible to find a function whose Laplace transform is any given rational function of s, by converting into a sum of partial fractions and then finding functions whose Laplace transform is each partial fraction.
The Inverse Laplace Transform
allows us to translate from the s domain (or frequency domain) back to
the t domain, or time domain. While there is a definition of the inverse
transform as a complex integral, in practice it is more useful to
derive it from the table of forward transforms, together with
the completion of the square and the following two properties of
the forward transform in which, as usual
.
{\cal L}(eatf(t)) = F(s-a)$$
This result is readily obtained by change of the variable of integration.
This result can be shown by differentiating under the integral sign.
With these results under our belt we can specify the following steps to finding the inverse transform for a rational function:-
Where N(s) and D(s) are polynomials in s. We may assume that the highest power of s in D(S) is of the form sn (for if not we can make it so by dividing N(s) and D(s) by the coefficient of sn in D(s)). \begin{enumerate} \item Factor D(s) into linear complex factors
In general this step will have to be performed numerically, and requires
D(s) to have numeric coefficients. Some algebraic manipulation programs
(eg. Macsyma) provide symbolic factorisation algorithms, but they cannot
always work (Galois). Note that finding eigenvectors of matrices also
requires the factorisation of polynomials, and it should therefore be no
surprise that a matrix formulation of the problem of solving Linear ODE's
with constant coefficients exists.
Thus D(s) will have factors of the form s-ai. If ai is complex,
then there will be a corresponding factor of the form s-\bar{ai},
where \bar{ai} is the complex conjugate of ai.
\item Combine the complex conjugate factors together to form real quadratic
factors. Thus
Where some of the factors may be repeated, and all of the a'i are real. In practice these first two steps will usually be combined. \item Express F(s) as partial fractions.
\item Form f(t), the inverse transform, as
f(t) = \sum fi(t) + \sum gi(t) $$
where the fi are derived by from the terms with linear denominators, the gi from those with quadratic denominators. Then, if ni = 1 :-
fi(t) = pi e^{a'it}$$ or if ni>1 :-
fi(t) = \frac{pi}{(ni-1)!} t^{ni-1} e^{a'it} $$
Proceeding to the quadratic terms we first complete the square in the denominator, obtaining
gi(t) = {\cal L}-1 \left(\frac{qis+ri}{(s2+bis+ci)^{mi}}\right)$$
= \invlap\left(\frac{qis+ri} {(s2+bis+(bi/2)2+(ci-bi2/4))^{mi}}\right)$$
= \invlap\left(\frac{qi(s+bi/2)+ri-qibi/2} {((s+bi/2)2 + di2)^{mi}}\right) = e^{-bit/2}\invlap\left(\frac{qis+r'i}{(s2+di2)^{mi}}\right)$$
Where di is real, since otherwise the denominator would factorise into real linear factors. The final step involves using the shifting theorem. Here ri' = ri-qibi/2
Now, if mi = 1:
gi(t) = e^{-bit/2}\invlap(\frac{q'is+r'i}{s2+di2})$$
= e^{-bit/2}(q'i\,\cos dit + \frac{r'i}{di}\sin dit) $$
If mi > 1 we can use the following results, which are given without justification:
\invlap\frac{1}{((s+a)2 + b2)n} = \frac{-e^{-at}}{4^{n-1}b^{2n}} \sum_{r=1}n \left(\begin{array}{c}2n-r-1\\n-1\end{array}\right) (-2t)^{r-1} \frac{d^r}{dt^r}\cos bt$$ and:
\invlap \frac{s}{((s+a)2 + b2)n} = \frac{e^{-at}}{4^{n-1}b^{2n}}\left\{ \sum_{r=1}n \left(\begin{array}{c}2n-r-1\\n-1\end{array}\right) \frac{1}{(r-1)!} (-2t)^{r-1}\frac{d^r}{dt^r}(a \cos bt + b \sin bt) \right.$$
\left. -2b \sum_{r=1}^{n-1} \left(\begin{array}{c}2n-r-2\\n-1\end{array}\right) \frac{1}{(r-1)!} (-2t)^{r-1}\frac{d^r}{dt^r}(\sin bt)\right\}$$ \end{enumerate}
Let us now consider a simple controller used this time to control a motor which actuates a joint of a 1-link robot. Recall that the thermostat was unable to control the room temperature very precisely. Perhaps if we used a continuously variable amplifier instead of a discrete switch we might get better control. This is indicated in the following diagram:-
Here the amplifier is supposed to be a current amplifier, ie one
whose output current is specified by the input (this simplifies our
description of the motor).
Now the motor produces a torque
related to the current i
by
where cM is a constant of the motor. Let the inertia of the link be J, and the angle of the joint be y. Then the motion of the joint is described by the angular analog of Newton's first law:-
Now the current i is commanded by the controller as follows:-
i = g(x-y)$$
Where g is the gain. (we assume that the amplifier converts 1 unit of output from the gain to 1 unit of current). Thus we have the differential equation for y
For convenience we let \omega2=\frac{gcM}{J}. This gives
\ddot{y} = \omega2(x-y)$$
If x=0 this will be recognised as the equation of simple harmonic motion - so it seems that our controller will oscillate, which is not a desirable state of affairs. Let us verify this by transforming the equation into the s-domain.
Solving for Y(s) we get:-
Thus if we choose x(t) to be the unit step u(t), we get
Expanding as partial fractions we get:-
From which we see, by equating coefficients, that:-
Transforming back into the t-domain, we obtain:-
y(t) = u(t) - \cos\omega t$$
That is to say, y differs from x by the term \cos\omega t, so that the link oscillates for ever about the goal position. We can see how to improve the controller by considering a pendulum, or mariner's compass, whose motion is damped by being placed in fluid. The fluid exerts a viscous friction which slows down the oscillations. For certain conditions it is reasonable to treat this damping as being proportional to the velocity. (You may know that sub-sonic aerodynamic form drag is taken to be proportional to the square of the velocity).
Likewise we can damp out the oscillations of our controlled link by subtracting off a term proportional to the velocity. The velocity can either be found from a tachomoter, or by differentiating the position (but this is a noisy estimate).
\begin{picture}(5,2)
\put(1,1){\framebox[0.5in]{\Sigma}}
\put(1.5,1.05){\vector(1,0){0.5}}
\put(2,1){\framebox[0.5in]{g}}
\put(2.5,1.05){\vector(1,0){0.5}}
........
\end{picture}
In the new improved controller the current i is commanded by the controller as follows:-
i = g(x-y) - k\dot{y}$$
Where k specifies the damping. Thus we now have differential equation for y
\ddot{y} = \frac{gcM}{J}(x-y) -\frac{kcM}{J}\dot{y} $$ As before \omega2=gcM/J, and we now choose \xi so that:-
2\omega\xi=\frac{kcM}{J}$$ This gives
\ddot{y} = \omega2(x-y)-2\omega\xi\dot{y}$$
Transforming into the s-domain as before:-
s2Y(s) = \omega2(X(s)-Y(s))-2\omega\xi sY(s)$$
Solving for Y(s) we get:-
Y(s) = \frac{\omega2} {s2+2\omega\xi s+\omega2}X(s)$$
Thus if we choose x(t) to be the unit step u(t), we get
Y(s) = \frac{\omega2}{(s2+2\omega\xi s+\omega2)s} $$
Expanding as partial fractions we get:-
Y(s) = \frac{A}{s} + \frac{Bs+C}{s2+2\omega\xi s+\omega2}$$
= \frac{(A+B)s2 + (C+2A\omega\xi)s + A\omega2} {s(s2+2\omega\xi s+\omega2)}$$
From which we see, by equating coefficients, that:-
A=1,\;\;\;\;C= -2\omega\xi,\;\;\;\;B=-1$$
Completing the square, the quadratic denominator becomes:-
s2+2\omega\xi s + (\omega\xi)2 + \omega2(1-\xi2)$$
Here, if \xi = 0 we have the undamped case we analysed earlier, if \xi<1 we have the underdamped case, if \xi=1 we have the critically damped case, and if \xi>1 we have the overdamped case.
In the underdamped case, let \omega'2=\omega2(1-\xi2). The behavior in the t-domain is then
y(t) = u(t) - \invlap\left( \frac{s + \omega\xi + \omega\xi}{(s+\omega\xi)2+\omega'2}\right)$$
And, applying the shifting theorem
y(t) = u(t) - e^{-\omega\xi t} \invlap\left(\frac{s+\omega\xi}{s2+\omega'2}\right)$$
= u(t) + e^{-\omega\xi t}(-\cos \omega't -\frac{\omega\xi}{\omega'}\sin\omega't)$$
In the critically damped case, we have a perfect square down below, so we can either apply the tn rule, or treat it as a limiting case of underdamped behavior, remembering that \cos(x)\rightarrow 1 and sin(x)/x \rightarrow 1 as x\rightarrow 1. Thus we obtain, for the critically damped case:- NOTE write it as a limiting condition.
y(t) = u(t) - e^{-\omega t} - \omega t e^{-\omega t}$$
Now as t\rightarrow \infty, t e^{-\omega t} \rightarrow 0 (exponentials increase faster than any polynomials.), so as t\rightarrow \infty, y(t)\rightarrow 0.
So y(0+) = 1 - 1 - 0 = 0 i.e. the system starts at y=0
y'(t) = \omega e^{-\omega t} - \omega e^{-\omega t} + \omega2 t e^{-\omega t} = \omega2 t e^{-omega t} $$
Thus y'(0+) = 0 , so the system starts at rest, and y'(t) >0 when t>0, so that y(t) has no maximum value, but increases monotonically to the limit 1 as t tends to infty, i.e. there is no overshoot.
In the overdamped case, the denominator will factor into two real linear terms in s, so that the behavior in the t-domain will be an exponential approach to u(t). The \sinh and \cosh laplace transforms provide the neatest derivation. Let \omega"2=\omega2(1-\xi2). The behavior in the t-domain is then
y(t) = u(t) - \invlap\left( \frac{s + \omega\xi + \omega\xi}{(s+\omega\xi)2-\omega"2}\right)$$
And, applying the shifting theorem
y(t) = u(t) - e^{-\omega\xi t} \invlap\left(\frac{s+\omega\xi}{s2-\omega"2}\right)$$
= u(t) + e^{-\omega\xi t}(-\cosh \omega"t -\frac{\omega\xi}{\omega"}\sinh\omega"t)$$
Now the \cosh(t)\rightarrow\infty, \sinh(t)\rightarrow\infty as t\rightarrow\infty, so we have to satisfy ourselves that the system is stable. Using the definition cosh(x) = (e^x+e^(-x))/2
\begin{Large}NOTE \end{Large} If two functions have the same Laplace transform they must have the same value ``almost everywhere" on the positive real axis. This means that they may differ at finite number of points (or more depending on the definition of integration employed, eg. is it Lebegue Integration)
The control systems we have studied so far have had the property that, provided they are stable, if they are given a step input, they approach the new commanded value asymptotically, ie if x(t) = au(t) then \lim_{t \rightarrow \infty} y(t) = 1 This behavior is clearly desirable for a robot, but it is not always attainable. Suppose for example that our link is loaded by gravity. Then in order for the link to be stationary, the motor must exert a torque to resist the gravitational force, and therefore there must be a non-zero current in the motor. But i = g(x-y)+k\dot{y} so that a non-zero current can only exist with a stationary link if there is a non-zero error. Now if g is large, the error will be small, but we have remarked that we may want to position a robot to one part in 10000, and that we cannot have too large a gain or the system will be unstable. (Also if we are digitising our sensor input, we will not be able to position our robot exactly to one quantisation level, since to get any current we will have to have at least one quantisation level of error). There is more than one solution to this problem. For example we could :- \begin{enumerate} \item Calculate what the current required to overcome gravity is, and add this into the command to the motors. In the case of gravity, this quantity depends only on the joint angle y, so call it p(y) Then we use the equation
i = g(x-y) +k\dot{y}+ p(y)$$ If we get the sums right, then gravitational loads no longer produce any error. The approach can be extended to handle dynamic loads arising from trying to track a changing input (eg a ramp x(t)=a\,r(t)). However it will not work for friction, since the value of frictional forces cannot be predicted with any accuracy, and there is an inconvenience if the robot is carrying a load, since the static (and dynamic) characteristics of the load have to be known. \item Consider the integral of the error with respect to time. If the error is of a single sign, and non-zero, then its integral is a monotonic increasing function. If x(t)-y(t)>\epsilon>0 where \epsilon is a constant, then \int(x-y)dt will increase without limit. So if we employ a controller for which i = gP(x-y) - gV\dot{y} + gI\int(x-y)dt, we can be sure that if we give it a step command, the error will not remain greater than some constant \epsilon. It may, of course, oscillate, and will do so if we choose too large a value for gI. \end{enumerate}
A controller which employs Position, its Derivative, and Integration is known as a PID controller. The earlier controller which employed just the position and derivative of y is known as a PD controller.
The limiting behavior of a linear controller can be derived using the following Final Value Theorem.
If F(s) = {\cal L}(f(t)) then
\lim_{t\rightarrow\infty} f(t) = \lim_{ s\rightarrow 0} sF(s)$$
provided that F(s) is analytic in the right half plane including the j axis. (For a rational F(s) this means that it must have no poles in the right half plane). Thus, applying this theorem to our PD controller with a unit step input, we get:-
sY(s) = \frac{\omega2}{(s2+2\omega\xi s+\omega2)} $$
So:-
\lim _{t \rightarrow\infty} y(t) = \lim _{s \rightarrow 0} sY(s) = 1 $$
If now we have a disturbing force f(t), then
\ddot{y} = \frac{gcM}{J}(x-y) -\frac{kcM}{J}\dot{y}+f/J $$
Choosing \omega and \xi as before, we have:-
\ddot{y} = \omega2(x-y)-\omega\xi\dot{y}+f/J$$
Transforming into the s-domain as before, with F(s) as the laplace transform of f(t):-
s2Y(s) = \omega2(X(s)-Y(s))-\omega\xi sY(s) + F(s)/J$$
Solving for Y(s) we get:-
Y(s) = \frac{\omega2X(s) + F(s)/J} {s2+2\omega\xi s+\omega2}$$
Thus if we choose x(t) to be the unit step u(t), and f(t) to be the unit step f1u(t) we get:-
Y(s) = \frac{\omega2/s +f1/(sJ) }{(s2+2\omega\xi s+\omega2)} $$
(For the purpose of the analysis, we have to assume that the disturbing force starts being applied at t=0 because of the requirement that all variables be zero initially).
Thus
\lim _{t\rightarrow\infty} y(t) = 1 + f1/(\omega2J)$$
So that we have a steady state error.
Analysing the system with an integral term, we obtain:-
s2Y(s) = \omega2(X(s)-Y(s)) + gI(X(s)-Y(s))/s-\omega\xi sY(s) + F(s)/J$$
s^3Y(s) = s\omega2(X(s)-Y(s)) + gI(X(s)-Y(s)) - \omega\xi s2Y(s) + sF(s)/J$$
Y(s) = \frac{(s\omega2+gI)X(s) + sF(s)/J}{s^3+s2\omega\xi+s\omega2+gI}$$
If x(t) = u(t) then ,applying the final value theorem, we see that y(t) tends to 1 as t\rightarrow\infty, provided the denominator has no roots in the right half plane. Now if gI is large, we can see thet the denominator behaves as s^3+gI which has roots s=\alpha e^{j(\pi/3+2k\pi/3)} where \alpha^3=gI. k=1 gives the root s=-\alpha, while k=0 and k=2 give two roots in the positive half plane. Thus for a large gI we may expect the system to be unstable. If gI = 0 we have a root s=0 and the two roots of the denominator of the system analysed earlier, which will be complex in the underdamped case. As gI is increased from zero, the real root will start to move along the negative real axis, and the complex roots will initially remain in the left half plane, so we have a stable system for small values of gI (to show that the real root moves this way, hold all coefficients constant except gI and consider \frac{ds}{dgI}, which is negative if s=0).
The position of robot joints is usually controlled digitally since this provides greater accuracy and flexibility. However, unless great computer power is used to give short sample times, it may be advantageous to provide velocity feedback through an analog controller. Recall that a long sample time can mean instabilility. If we restructure the equation
i = g(x-y) + k\dot{y}$$ as
v= g1(x-y),\;\;\;\; i = k(v-\dot{y})$$
Where g1 = g/k We can regard the second equation as defining a velocity controller which is fed a velocity command by the first. This arrangement is more stable if the velocity control is performed by an analog controller which by its nature imposes no significant delay, and the position control is performed by a digital controller. The precision of velocity control does not have to be great for most robotic use, since precision of position is normally what is required, so that the analog velocity loop is very adequate from that point of view.
While digital controllers can readily represent a very large range of numbers, the analog components of a system, the power amplifiers and the motors, will only have a limited range in which their behavior is linear. If the input to an amplifier goes outside the operational range the device will saturate, that is to say it will give a constant output. A motor (or at least vital parts of it) may well melt. In practice, the gains of a controller are such, that if a step input is given to move any significant fraction of the full range of movement of a joint, the amplifiers will saturate. This may be acceptable in a single link device, provided the amplifiers are properly designed to go into saturation, and properly matched to the motors. However once a controller goes into saturation, it is no longer controlling. And if we have multiple joints on a robot, they will behave in a very uncoordinated manner if we give them all large step inputs (eg. to move the robot from one place to another), so that the path taken by the device will be unpredictable.
A practical solution to this problem, suitable for moderate performance robots, is to interpolate in joint space. If we have joints with angular or linear displacements (y1,y2,...,yn) then we can choose command functions (x1(t)...xn(t)) to go from a point (a1...an) to (b1...bn) by a straight line.
\vec{x(t)} = (1-\lambda(t))\vec{a} + \lambda(t)\vec{b}$$
There remains a problem of choosing the function \lambda. One possibility is to have \lambda(t)=lt for some constant l, which can be chosen by the programmer of the robot. He will typically develop the program with a small value of l and then increase it until the robot starts behaving improperly. This is not perhaps a very scientific approach, but it is often done. If we have a kinematic transform from cartesian space to joint space, the contoller can work out a straight line in cartesian space, and transform that to a sequence of values for x(t). Note that there is no guarantee that y(t) will track x(t) exactly. Indeed, if \lambda(t)=lt then it is impossible to obtain exact tracking, since a robot cannot start moving instantaneously. However this form of interpolation does provide a practical approach to controlling a robot which giving step inputs does not.
One of the most important aspects for the controllability of a robot is what the lowest resonant frequency of vibration is. One can think of a robot as a device that converts a control signal to a movement. Consider for example a robot that is to pick up an object at point A and put it down at point B, and repeat the action every second. If it has a resonant frequency below 2hz it will be very difficult to control it to perform this task.