Differential flatness for neuroscience population dynamics -- A preliminary study
HHugues Mounier
Differential flatness for
Neurosciencepopulationdynamics
A preliminary study
Version – November a r X i v : . [ c s . S Y ] D ec i ii Hugues MounierLaboratoire des Signaux et SystèmesCentraleSupélec , rue Joliot Curie GIF sur YVETTEe-mail: [email protected]
O N T E N T SI neural mass models 31 simple and common neural mass models 51 . Scalar integrate and fire models . . . . . . . . . . . . . . . . . . Two variables integrate and fire models . . . . . . . . . . . . . Two variables integrate and fire Izhikevich’s models . . . . . . Vectorial integrate and fire models . . . . . . . . . . . . . . . . Neural mass Wilson-Cowan E-I networks . . . . . . . . . . .
72 differential flatness 92 . Differential flatness notion . . . . . . . . . . . . . . . . . . . . . . Dynamics and observation equations . . . . . . . . . . . Differential flatness definition . . . . . . . . . . . . . . . . Parametrization . . . . . . . . . . . . . . . . . . . . . . . . A word of methodology . . . . . . . . . . . . . . . . . . Differential flatness of simple neural mass models . . . . . . . . Weakly coupled E-I networks . . . . . . . . . . . . . . . . Differential flatness of a weakly coupled E-I network . . Differential flatness of Wilson Cowan’s E-I network . . . Differential flatness of asymetric Wilson Cowan’s E-Inetwork . . . . . . . . . . . . . . . . . . . . . . . . . .
123 diff . flatness applications & extensions 153 . Differential flatness applications . . . . . . . . . . . . . . . . . . . Trajectory tracking . . . . . . . . . . . . . . . . . . . . . . Feedforward to feedback switching . . . . . . . . . . . . Cyclic character . . . . . . . . . . . . . . . . . . . . . . . . Positivity & Boundedness . . . . . . . . . . . . . . . . . . Simultaneous synchronisation & tracking . . . . . . . . Flatness appls. for neural mass E-I networks . . . . . . . . . . . . Trajectory tracking for weakly coupled E-I networks . . Trajectory tracking for asymetric Wilson-Cowan’s E-Inetworks . . . . . . . . . . . . . . . . . . . . . . . . . . . . Feedforward to feedback switching . . . . . . . . . . . Diff. flatness appls. for simplistic motor control . . . . . . . . . . Two link arm model . . . . . . . . . . . . . . . . . . . . . Differential flatness and open loop control of the twolink arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . End effector dynamics . . . . . . . . . . . . . . . . . . . . Trajectory tracking of the two link arm . . . . . . . . . Extensions of differential flatness . . . . . . . . . . . . . . . . vi Contents . Flatness of other neural mass models . . . . . . . . . . . . . . . . Jansen and Rit model . . . . . . . . . . . . . . . . . . . II neural field population models 373 . . General case model . . . . . . . . . . . . . . . . . . . . . . Parameters and variables assumptions . . . . . . . . . . . Neuron interaction strength examples . . . . . . . . . . . Simplification hypotheses . . . . . . . . . . . . . . . . . . Pointwise neuron interaction strength models . . . . . . Exponential type neuron interaction strength . . . . . . A Jansen and Rit Neural field model . . . . . . . . . . . . . .
III appendix 51a notations , transforms and sigmoid functions 53a . Functions and distributions . . . . . . . . . . . . . . . . . . . . . Transforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Sigmoid functions . . . . . . . . . . . . . . . . . . . . . . . . .
54b some flatness simple criteria 57b . Necessary and sufficient conditions in peculiar cases . . . . . . A necessary condition . . . . . . . . . . . . . . . . . . . . . . . . Static state feedback linearizability criterion . . . . . . . . . . . . Brief recall of differential geometry notions . . . . . .
58c definitions for extensions of differential flatness 61c . Systems, dynamics and differential flatness . . . . . . . . . . . . Basic definitions from differential algebra . . . . . . . . . Algebraic and transcendental extensions . . . . . . . . . Nonlinear systems and flatness . . . . . . . . . . . . . . H-fields, Liouvillian and existential closedness . . . . . . . .
62d two link arm inverse kinematics 65e two link arm code listing 67e . Flatness based control of the arm . . . . . . . . . . . . . . . . ontents notations Some recurrent notations, transforms and sigmoid functions we shallbe using throughout this document are given in Appendix A p. . introduction The various objectives one whishes to attain through a controlled dy-namical system almost always boil down to a system’s behavior mod-ification. One has at his disposal so called control variables whoseaim is to steer the system. The behavior modification generally con-sist in tracking a prescribed trajectory, with stability. Note that thefirst phase (open loop trajectory tracking) is a feedforward one, whilethe second (with stability) is a feedback one. A possible methodologyfor controlling a system is then decomposed in two steps: . A so called open loop trajectory tracking, supposing the modelperfect and the initial conditions perfectly known. . A feedback stabilizing the system around the reference trajecto-ries, to compensate for model mismatch, poorly known initialconditions and external perturbations.Simple and natural solutions to problem are obtained through thedifferential flatness property, a notion due to Michel Fliess, Jean Lévi-ne, Philippe Martin and Pierre Rouchon (Fliess et al., ). This prop-erty amounts to a parametrization of a dynamical system in terms ofa so-called flat output: any variable in the system can be expressedthrough a function of the flat output components and a finite numberof its derivatives. This parametrization yields expressions of all thesystem’s variables without having to integrate any differential equa-tion, which ensures fast computations. The control is in particulargiven through the flat output and its derivatives ; an inversion of thesystem control input to flat output is thus performed, without anyintegration. This notion has been extended to infinite dimensionalsystems, governed by partial differential equations (Woittennek andMounier, ).The present document is devoted to structural properties of neuralpopulation dynamics and especially their differential flatness. Severalapplications of differential flatness in the present context can be en-visioned, among which: trajectory tracking, feedforward to feedbackswitching, cyclic character, positivity and boundedness. art IN E U R A L M A S S M O D E L S S I M P L E A N D C O M M O N N E U R A L M A S S M O D E L SThe following models are quite simple models for neural populationswith lumped space parameters (see, e.g. (Dayan and Abbott, ),Chapter , (Ermentrout and Terman, ), Chapter ). The case ofdistributed parameter models, so-called neural field models, is con-sidered in Section . . , p. . . These type of models are of the form: C ˙ ν = − g L ( ν − ν L ) + F ( ν ) + I ( . )where ν is the membrane potential, determined with respect to theresting potential of the cell, τ m is the membrane time constant, F ( ν ) is a spike generating current, and I is the total current elicited bysynaptic inputs to the neuron.Common types of models, each associated with a specific type ofspike generating currents, are: ◦ The leaky integrate and fire, corresponding to F = ◦ The quadratic integrate and fire (or theta neuron), correspond-ing to F ( ν ) = g L ∆ T ( ν − ν T ) + g L ( ν − ν T ) − I T ◦ The exponential integrate and fire, corresponding to F ( ν ) = g L ∆ T e ν − ν T ∆ T . More general models adds a second variable coupled to the voltage(see (Izhikevich, ))˙ ν = F ( ν ) − µ + I ˙ µ = a ( b ν − µ ) The function F ( ν ) describes the current–voltage characteristic of themembrane potential near the threshold, and it typically looks like aparabola (Izhikevich, , ): F ( ν ) = ν . Other choices possibleare F ( ν ) = | ν | , F ( ν ) = − ν , F ( ν ) = | ν | n + − ν An exponential spike generating current has been considered in (Bretteand Gerstner, ) leading to the so-called adaptive exponential in-tegrate and fire: F ( ν ) = e ν − ν , and (Touboul, ) suggested thequartic model F ( ν ) = ν + a ν . . ’ s mod - els Another class of models is found in (Izhikevich, ):˙ ν = F ( ν ) − µ ( E − ν ) + I ˙ µ = a ( b ν − µ ) where ν plays the role of a conductance and E is its reverse potential,which could be assumed to take values ± . Two types of integrate and fire models can be derived (see, e.g. (Er-mentrout and Terman, ), Chapter ): τ m ˙ ν = − ν + W F ( ν ) + ˜ I ( . )where τ m is the membrane time constant, and τ d ˙ ρ = − ρ + F ( W ρ + I ) ( . )where τ d is the synaptic decay time. . - cowan e - i networks 7 Remarks . The choice of one of the models is based on time scale con-siderations (see, e.g. (Ermentrout and Terman, ), p. ), wherein ( . ) the temporal dynamics is dominated by the synaptic decay andin ( . ) , the membrane time constant of the postsynaptic cell are smallcompared with the decay of the synapse. . Note that the above two models can be shown to be equivalent (when τ d = τ m = τ ) in the following sense (see (Miller and Fumarola, ).If ρ is a solution of the membrane model ( . ) , then W ρ + I is asolution of ( . ) . Indeed, setting ν = W ρ + I , one obtains τ ˙ ν = τ W ˙ ρ + τ ˙ I = W (cid:0) − ρ + F ( W ρ + I ) (cid:1) + τ ˙ I = − ( ν − I ) + W F ( ν ) + τ ˙ I = − ν + W F ( ν ) + ˜ I . - cowan e - i networks Consider the simplest form of network, a pair of mutually coupledlocal populations of excitatory and inhibitory neurons, also called E-Inetwork (see, e.g. (Bressloff, ), Subsection . , p. ). This modelwas originally developed by Wilson and Cowan (see, e.g. (Ermentroutand Terman, ), Subsection . , p. ), and has the form τ e ˙ ν e = − ν e + F e ( w ee ν e − w ie ν i + I e ) ( . a) τ i ˙ ν i = − ν i + F i ( w ii ν i − w ei ν e + I i ) ( . b)where ν e and ν i are the proportion of excitatory and inhibitory cellsfiring per unit time, the activations are nonlinear functions (typicallysigmoidal) F e , F i of the presently active proportion of cells, w ∗ are thestrength of the connections.The matrix form of the previous model is τ e ˙ ν e = − ν e + F e ( W ee ν e − W ie ν i + I e ) ( . a) τ i ˙ ν i = − ν i + F i ( W ii ν i − W ei ν e + I i ) ( . b) D I F F E R E N T I A L F L AT N E S S . . . Dynamics and observation equations
Consider a system given by the dynamics equation and the observa-tion equation˙ x = f ( x , u ) dynamics equation ( . a) y m = h ( x ) observation equation ( . b)with x ( t ) = ( x ( t ) , . . . , x n ( t )) , the state , or, in Karl Friston’s terms the hidden variables (see, e.g. (Friston, )), i.e. the controlled variables, u ( t ) = ( u ( t ) , . . . , u m ( t )) , the control input , functions enabling an ac-tion on the process (typically input current), and y m = ( y m ( t ) , . . . , y mp ( t )) the output , measured functions enabling to sense the environment(quantities coming from sensors).Note that the dynamics equations form an undetermined system ofdifferential equations, since the control functions u ( t ) are not a pri-ori determined. Once the control variables are fixed (i.e. substitutedwith known functions of time), the system ( . ) becomes determined(i.e. can be integrated). The state variables represent the instanta-neous memory of the system: once the control variables have beendetermined, the knowledge of the state variables (at time t ) enablesto predict the future state (at time t + dt ).Another formualtion is the following: the state of a dynamical sys-tem is a set of physical quantities the specification of which (in theabsence of external excitation) completely determines the evolutionof the system. . . Differential flatness definition
The notion of differential flatness (see (Fliess et al., )) is a formof controllability for non linear dynamical systems which is espe-
90 differential flatness cially well suited for trajcetory tracking problems. It amounts to aparametrization of the system without integration of any differentialequation. Although the mathematical property seems quite strong, itappears that this notion is commonly encountered in practice (see, e.g.(Rouchon, , Martin and Rouchon, ) for a catalog of differen-tially flat systems). We shall give below a definition for such systemsand illustrate this through simple examples derived from the wellknown Wilson and Cowan’s model. Some more details about thisproperty is given in the appendices.
Definition The system ˙ x = f ( x , u ) ( . ) with x ( t ) ∈ R n and u ( t ) ∈ R m is differentially flat if there exists a set ofvariables, called a flat output , y = h ( x , u , ˙ u , . . . , u ( r ) ) , y ( t ) ∈ R m , r ∈ N ( . ) such that x = A ( y , ˙ y , . . . , y ( ρ x ) ) ( . a) u = B ( y , ˙ y , . . . , y ( ρ u ) ) ( . b) with q an integer, and such that the system equationsdAdt ( y , . . . , y ( q + ) ) = f ( A ( y , . . . , y ( q ) ) , B ( y , . . . , y ( q + ) )) are identically satisfied. . . Parametrization
For any flat output given through a function of the form t ∈ R → y ( t ) ,the trajectory of the system x ( t ) , u ( t ) are given by: x ( t ) = A ( y ( t ) , ˙ y ( t ) , . . . , y ( ρ x ) ( t )) ( . a) u ( t ) = B ( y ( t ) , ˙ y ( t ) , . . . , y ( ρ u ) ( t )) ( . b)There is a one to one correspondance between the system trajectoriesand the ones given by the flat output. . . . A word of methodology
The preceding notion will be used to obtain so called “open loop”controls, that is control laws which will ensure the tracking of the ref-erence flat outputs when the model is assumed to be perfect and the stateinitial conditions are assumed to be exactly known . Since this is never thecase in practice, one needs some feedback schemes that will ensureasymptotic convergence to zero of the tracking errors. Our frame-work can thus be decomposed in two steps: . Design of the reference trajectory of the flat outputs; off-linecomputation of the open loop controls. . Inline computation of the complementary closed loop controlsin order to stabilize the system around the reference trajectories.Why is this two step design better suited than a classical stabi-lization scheme? The first step obtains a first order solution to thetracking problem, while following the model instead of forcing it (likein a usual pure stabilization scheme). The second step is a refine-ment one, and the error between the actual values and the trackedreferences will be much smaller than in the pure stabilization case. . . . Weakly coupled E-I networks
Consider a Wilson-Cowan model where w ie (cid:28)
1. Hence ( . a) re-duces to τ e ˙ ν e = − ν e + F e ( w ee ν e + I e ) Alternatively, one could also consider the other limit case where w ei (cid:28) . b) reduces to τ i ˙ ν i = − ν i + F i ( w ii ν i + I i ) Whatever case we consider, we shall abbreviate it by the followinghighly simplified model τ ˙ ν = − ν + F ( w ν + I ) ( . ) where the subscript has been dropped for convenience. This model,although simplistic, is considered here because of its simplicity forpedagogical purposes. Set φ = F − where F is a sigmoid function (see A. , p. ). . . Differential flatness of a weakly coupled E-I network
The model depicted by ( . ) is differentially flat, with ν as a flat out-put. Indeed, one has w ν + I = φ ( τ ˙ ν + ν ) and the input I is given by I = − w ν + φ ( τ ˙ ν + ν ) ( . ) . . Differential flatness of Wilson Cowan’s E-I network
The E-I network equations ( . a)–( . b) τ e ˙ ν e = − ν e + F e ( w ee ν e − w ie ν i + I e ) τ i ˙ ν i = − ν i + F i ( w ii ν i − w ei ν e + I i ) rewrite w ee ν e − w ie ν i + I e = F − e ( τ e ˙ ν e + ν e ) w ii ν i − w ei ν e + I i = F − i ( τ i ˙ ν i + ν i ) This model is thus differentially flat with flat output ( ν e , ν i ) : I e = w ie ν i − w ee ν e + F − e ( τ e ˙ ν e + ν e ) I i = − w ei ν e + w ii ν i + F − i ( τ i ˙ ν i + ν i ) . . Differential flatness of asymetric Wilson Cowan’s E-I network
Consider the E-I Wilson-Cowan equations ( . a)–( . b) with staticallycoupled external currents I e = ( + a ) I , I i = ( − a ) I . with a ∈ [ −
1, 1 ] an asymetry factor. The model ( . a)–( . b) then be-comes (see, e.g. (Ermentrout and Terman, ), Section . , p. ) τ e ˙ ν e = − ν e + F e (( + a ) I − w ie ν i ) τ i ˙ ν i = − ν i + F i (( − a ) I − w ei ν e ) In the asymetric case, i.e. a = − τ e ˙ ν e = − ν e + F e ( − w i ν i ) ( . a) τ i ˙ ν i = − ν i + F i ( I − w e ν e ) ( . b)Then, ν e is a flat output. Indeed, one gets ν i = − w i F − e ( τ e ˙ ν e + ν e ) I = (cid:104) w e ν e + F − i ( τ i ˙ ν i + ν i ) (cid:105) D I F F E R E N T I A L F L AT N E S S A P P L I C AT I O N S A N DE X T E N S I O N S . A number of applications of differential flatness can be envisioned,among which: ◦ Trajectory tracking . ◦ Feedforward to feedback switching . ◦ Cyclic character . ◦ Positivity & boundedness . ◦ Simultaneous synchronisation & tracking . . . Trajectory trackingFlatness and feedback linearization
A characterization of flat systems that appears very useful for stabi-lized trajectory tracking is the following
Proposition A system is flat if, and only if, it is linearizable by endoge-nous feedback and change of coordinates.
A dynamic feedback is called endogenous if it does not include anyexternal dynamics. More precisely
Definition Consider the dynamics ˙ x = f ( x , u ) . The feedback u = ξ ( x , z , v ) ˙ z = ζ ( x , z , v )
156 diff . flatness applications & extensions (where v is the new input) is called a dynamic endogenous feedback ifthe original dynamcis ˙ x = f ( x , u ) is equivalent to the transformed one ˙ x = f ( x , ξ ( x , z , v )) ˙ z = ζ ( x , z , v ) Two systems are called equivalent if there exists a invertible transformationwich exchanges their trajectories.
A more restrictive notion is the one of static state feedback, as de-scribed below.
Definition Consider the dynamics ˙ x = f ( x , u ) . The feedback u = ξ ( x , v ) (where v is the new input) is called a static feedback if the original dynam-cis ˙ x = f ( x , u ) is transformed to ˙ x = f ( x , ξ ( x , v )) See the Subsection B. , p. for a static state feedback linearizationcriterion. Dynamical extension algorithm
This procedure enables one to know if an m-uple ( y , . . . , y m ) is a flatoutput or not. Meanwhile, we shall obtain a linearizing feedback. phase i – Gathering the so called weak brunovsky indices. ) Differentiate y until a combination of controls appears. Note κ the number of successive differentiations y ( κ ) = f ) Differentiate y until a combination of controls (independent ofthe previous ones) appears. Note κ the number of successivedifferentiations y ( κ ) = f ...m) Differentiate y m until a combination of controls (independent ofthe previous ones) appears. Note κ m the number of successivedifferentiations y ( κ m ) m = f m . phase ii – Deciding the flatness character.Then, if κ + · · · + κ m = n ( n being the state dimension), the systemadmits ( y , . . . , y m ) as a flat output. If not, ( y , . . . , y m ) isn’t a flatoutput. phase iii – Obtaining the linearizing feedback.The linearizing feedback is given by f = v , . . . , f m = v m . Closed loop trajectory tracking
The open loop control laws suppose that the model is perfect andthat the initial conditions are exactly known. Since this is never thecase in practice, we add corrective terms to the open loop controlsderived above in order to stabilize the system around the referencetrajectories.More precisely, considering a flat dynamics ˙ x = f ( x , u ) , we wantto derive a controller able to follow any reference trajectory t (cid:55)→ y r ( t ) .In order to compensate for model mismatch and poorly known initialconditions, one has to complement the open loop (obtained throughflatness) with a closed loop corrective term depending on the error y ( t ) − y r ( t ) .Knowing the dynamics is flat, with flat output y , it can be trans-formed via endogenous feedback and coordinate change to a lineardynamics of the form y ( κ ) = v ... y ( κ m ) m = v m with the new input ( v , . . . , v m ) . Then, the elementary tracking feed-back v i = y ( κ i ) ir − κ i − ∑ j = k ij ( y ( j ) i − y ( j ) ir ) , i =
1, . . . , m = y ( κ i ) ir − κ i − ∑ j = k ij e ( j ) i with appropriately chosen k ij gains renders the error dynamics asym-totically stable: e ( κ i ) i = κ i − ∑ j = − k ij e ( j ) i , i =
1, . . . , m . flatness applications & extensions . . Feedforward to feedback switchingOpen and closed loop
The so-called open loop control u o is obtained through ( . b) u ( t ) = B ( y ( t ) , ˙ y ( t ) , . . . , y ( ρ u ) ( t )) by replacing y with a sufficiently differentiable trajectory y r ( t ) : u o ( t ) = B ( y r ( t ) , ˙ y r ( t ) , . . . , y ( ρ u ) r ( t )) The use of this control law would lead to the desired tracking behav-ior y = y r if the model ( . a) was perfect and if the initial conditionson y was precisely known. Since this is never the case in practice,one has to use closed loop feedback laws, such as the ones elaboratedin the previous Subsection . . , p. . The difference between openand closed loop control laws can be bounded by the tracking errorand its derivatives. The simple weakly coupled E-I network exampleis examined in Subsection . . , p. . Temporal switching from feedforward to feedback
Consider the following control law u ( t ) = ( − σ ( t − t sw )) u o ( t ) + σ ( t − t sw ) u c ( t ) ( . )with σ a sigmoid function, for example of the form σ ( t ) = + e − t − βα , σ ( t ) = + tanh ( α t ) t = t = t sw − d for some d >
0, we have u ≈ u o , andfrom t = t sw + d , u ≈ u c . This is the kind of control human beingstend to adopt for example in gesture control. When grasping a glassof water, the first part of the gesture is done in open loop, quicklyand inaccurately; the second part of it is done with visual feedback,much more slowly but precisely. . . Cyclic character
When the flat output is cyclic, i.e. ( − ∆ τ ) y ( t ) = y ( t ) − y ( t − τ ) = . Then, all the variables wich are expressed as functions of y ( t ) , that isall the variables when the system is flat, are also cyclic: for a variable z which is expressed as z = C ( y , ˙ y , . . . , y ( η ) )( − ∆ τ ) z = ( − ∆ τ ) C ( y , ˙ y , . . . , y ( η ) )= C (( − ∆ τ ) y , ( − ∆ τ ) ˙ y , . . . , ( − ∆ τ ) y ( η ) )= . )More generally, if the flat output satifies a difference equation: p ( ∆ τ ) y = p is a polynomial, then any variable of the flat system with flatoutput y also satifies the same difference equation. . . Positivity & Boundedness
The goal is here to specify the reference trajectory in order to enforcecertain properties for various system variables. Two main cases canbe considered: The one of cyclic reference trajectories and the one ofnon cyclic ones.
Cyclic reference trajectories
The flat output reference trajectories being cyclic can be expressedthrough a Fourier series ∀ i =
1, . . . , m , y i = ∞ ∑ n = ξ i , n e j π n ( . )Since all variables are also cyclic (see ( . )) they can also be expressedthrough a Fourier series z = ∞ ∑ n = ζ n e j π n ( . )One then has some relations expressing the ζ n through the ξ i , n : ζ n = φ n ( ξ n , . . . , ξ m , n ) And the positivity can be expressed through a sum of squares typeformula (see, e.g. (Dumitrescu, )). Several matlab packages areavailable for finding sum of squares decompositions of real multivari-ate polynomials; the most popular ones are SOSTOOLS (see http:// . flatness applications & extensions ), YALMIP (see http://users.isy.liu.se/johanl/yalmip/ , and especially http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Examples.MoreSOS ) and GloptiPoly (see http://homepages.laas.fr/henrion/software/gloptipoly/ ). Non cyclic reference trajectories
One option is to take in the flat output y = ( y , . . . , y m ) all the com-ponents y i s as polynomial splines. If the firing rate function is takento be of Naka Rushton type, then all inequalities will boil down toexpressions of the form P i ( y , ˙ y , . . . , y ( ρ ) ) > i =
1, . . . , m where the P i s are polynomials in their variables. Since the y i s arepolynomial splines, P i ( y , . . . , y ( ρ ) ) will be another polynomial spline.And any approximating polynomial spline is contained in the convexhull of its control points. One then can choose the lowest of these tobe positive, to ensure the above inequality to be fullfilled. . . Simultaneous synchronisation & tracking
One considers here two (or more generally N ) oscillators coupled viatheir input:˙ x = f ( x , u ) ˙ x = f ( x , u ) The flatness of this system ensures not only that synchronisation ispossible, but also that any periodic trajectory (of the flat output) maybe tracked, which is a much stronger result. . - ral mass e - i networks3 . . Trajectory tracking for weakly coupled E-I networks
Recall the weakly coupled E-I network ( . ), p. : τ ˙ ν = − ν + F ( w ν + I ) ( . ) . . for neural mass e - i networks 21 In a trajectory tracking, one chooses a reference trajcetory ν r andone wants that lim t → ∞ ν = ν r , or, what is the samelim t → ∞ e ν =
0, where e ν = ν − ν r This behavior can be enforced through the following desired errordynamics˙ e ν = − λ e ν ( . )where λ > . ), one has toset in ( . ): − ν + F ( w ν + I c ) = − τλ e ν + τ ˙ ν r where I c is the closed loop control law. The preceding equation canbe rewritten as w ν + I c = φ ( ν − τλ e ν + τ ˙ ν r ) which yields the following closed loop tracking feedback law I c = − w ν + φ ( τ ˙ ν r + ν − τλ e ν ) ( . )which ensures, through ( . ), the tracking of the reference trajectory ν r for the system ( . ) with stability. Remark The application of the preceding extension algorithm is quitetrivial since the system is fairly simple: ◦ Gathering the so called weak brunovsky indices.The flat output ν is differentiated once in equation ( . ) where thecontrol I is already present, hence κ = . ◦ Deciding the flatness character.Since the dimension of the state is n = , ∑ i κ i = κ = n and thesystem is flat with flat output ν . ◦ Obtaining the linearizing feedback.The linearizing feedback is given by: τ ( − ν + F ( w ν + I )) = v ( . ) . flatness applications & extensions This feedback transforms the dynamics ( . ) into the following linear one: ˙ ν = vand the elementary tracking feedback isv = ˙ ν r − λ e ν ( . ) Thus, the original tracking feedback law is obtained from ( . ) and ( . ) :I = − w ν + φ ( τ ˙ ν r + ν − τλ e ν ) . . Trajectory tracking for asymetric Wilson-Cowan’s E-I networks
Recalling the equations of the asymetric Wilson-Cowan’s E-I network( . a)–( . b) τ e ˙ ν e = − ν e + F e ( − w i ν i ) τ i ˙ ν i = − ν i + F i ( I − w e ν e ) and differentiating the first equation in ν e , we get the flat output dy-namics τ e ¨ ν e = − ˙ ν e − w i F (cid:48) e ( − w i ν i ) ˙ ν i = − ˙ ν e + w i τ i F (cid:48) e ( − w i ν i ) (cid:0) ν i − F i ( I − w e ν e ) (cid:1) ( . )The desired dynamics being¨ e er = − λ e e er − µ e ˙ e er , where e er = ν e − ν er the right hand side of ( . ) is then taken to be˙ ν e + w i τ i F (cid:48) e ( − w i ν i ) (cid:0) − ν i + F i ( I − w e ν e ) (cid:1) = τ e (cid:0) ¨ ν er + λ e e er + µ e ˙ e er (cid:1) Thus we get − ν i + F i ( I − w e ν e ) = τ i τ e w i F (cid:48) e ( − w i ν i ) (cid:0) τ e ˙ ν e + ¨ ν er + λ e e er + µ e ˙ e er (cid:1) and the tracking control feedback loop is obtained as I = (cid:104) w e ν e + F − i (cid:16) ν i + τ i τ e w i F (cid:48) e ( − w i ν i ) (cid:16) τ e ˙ ν e + ¨ ν er + λ e e er + µ e ˙ e er (cid:17)(cid:17)(cid:105) . . for neural mass e - i networks 23 . . Feedforward to feedback switchingOpen and closed loop
The so-called open loop control law is obtained when replacing ν bythe reference trajectory ν r in ( . ): I o = − w ν r + φ ( τ ˙ ν r + ν r ) ( . )The use of this control law would lead to the desired tracking behav-ior ν = ν r if the model ( . ) was perfect and if the initial conditionson ν was precisely known.This type of law is typically used by the brain, after training, forquick movements where the sensory system is bypassed. When thesensory system is used, the so-called closed loop control law ( . ) isapplied.The difference between I c and I o is I o − I c = w e ν − φ ( τ ˙ ν r + ν − τλ e ν ) + φ ( τ ˙ ν r + ν r )= w e ν − φ (cid:0) τ ˙ ν r + ν r + ( − τλ ) e ν (cid:1) + φ ( τ ˙ ν r + ν r ) Then, supposing φ to be globally γ -lipschitz: (cid:12)(cid:12) φ (cid:0) τ ˙ ν r + ν r +( − τλ ) e ν (cid:1) − φ ( τ ˙ ν r + ν r ) (cid:12)(cid:12) (cid:54) γ | − τλ || e ν | Hence the difference I c − I o admits the following bound | I c − I o | (cid:54) (cid:0) α + γ | − τλ | (cid:1) | e ν | Thus, if the tracking error is small, I c is in a neighborhhod of I o . Temporal switching from feedforward to feedback
Consider the following control law I ( t ) = (cid:0) − σ ( t − t sw ) (cid:1) I o ( t ) + σ ( t − t sw ) I c ( t ) ( . )with σ a sigmoid function (see A. , p. ). Thus, from t = t = t sw − d for some d >
0, we have I ≈ I o , and from t = t sw + d , I ≈ I c . This is the kind of control human beings tend to adopt forexample in gesture control. When grasping a glass of water, the firstpart of the gesture is done in open loop, quickly and inaccurately; the . flatness applications & extensions second part of it is done with visual feedback, much more slowly butprecisely. The expression can alternatively be expressed as: I = (cid:0) − σ sw (cid:1) φ ( τ ˙ ν r + ν r ) + σ sw (cid:0) − we ν + φ ( τ ˙ ν r + ν − τλ e ν ) (cid:1) = φ ( τ ˙ ν r + ν r ) + σ sw (cid:0) − we ν + ∆ φ ( ν , ν r ) (cid:1) where σ sw = σ ( t − t sw ) and ∆ φ ( ν , ν r ) = φ ( τ ˙ ν r + ν − τλ e ν ) − φ ( τ ˙ ν r + ν r ) . - tor control3 . . Two link arm model
Consider a two link robot arm acting as a simplistic model of a humanarm: M θ + M θ + C ( θ , ˙ θ ) + G ( θ ) = T ( . a) M θ + M θ + C ( θ , ˙ θ ) + G ( θ ) = T ( . b)where θ is the angle of the first arm, θ of the second, θ = ( θ , θ ) T , M ij are equivalent masses, C i are the coriolis forces, G i are the gravityforces, and T i are the control torques. The expressions for the C i ’sand the G i ’s are the following: ◦ The inertia expressions are M = J + J + m r + m ( l + r + l r cos θ ) ( . a) M = M = J + m ( r + l r cos θ ) ( . b) M = J + m r ( . c)where J i is the inertia of link i , m i its mass, l i its length, and r i thedistance from the beginning of the link to its center of mass. ◦ The Coriolis terms are given by: C = − m l θ r sin θ − m l ˙ θ ˙ θ r sin θ ( . a) C = m l ˙ θ r sin θ ( . b) ◦ And the gravity terms are G = ( m l + m r ) g sin θ + m r g sin ( θ + θ ) ( . a) G = m r g sin ( θ + θ ) ( . b) . . flatness appls . for simplistic motor control 25 Equations ( . ) can be rewritten in a vectorial form; to this purpose,set M = (cid:32) M M M M (cid:33) , θ = (cid:32) θ θ (cid:33) C = (cid:32) C C (cid:33) , G = (cid:32) G G (cid:33) , T = (cid:32) T T (cid:33) Then, model ( . ) becomes M ¨ θ + C ( θ , ˙ θ ) + G ( θ ) = T ( . ) θ l θ l Figure . : A two link robot arm. . . Differential flatness and open loop control of the two link arm
Remark Model ( . ) is differentially flat, with θ as a flat output. Indeed,the inputs T are directly expressed in terms of θ and its derivatives: T = M ¨ θ + C ( θ , ˙ θ ) + G ( θ ) and the open loop control for a trajectory θ r , θ r given by T r = M ¨ θ r + C ( θ r , ˙ θ r ) + G ( θ r ) Knowing that the desired trajectory is generally not given in termsof θ , θ but in terms of the end effector coordinates h x , h y , we have toexpress the former in terms of the latter. The end effector (e.g. thewrist) coordinates are given by: h x = l cos θ + l cos ( θ + θ ) ( . a) h y = l sin θ + l sin ( θ + θ ) ( . b) . flatness applications & extensions The inversion of these formulae is detailed in Appendix D, p. . Weshall here give the final expressions: θ = arctan (cid:18) h y h x (cid:19) − arctan (cid:18) l sin θ l + l cos θ (cid:19) ( . a) θ = arctan (cid:32) ± √ − ¯ h ¯ h (cid:33) ( . b)¯ h = h x + h y − l − l l l . . End effector dynamics
The dynamcis in θ is given by:¨ θ = − M − ( C + G ) + M − T ( . )and the dynamcis in the end effector, i.e. in h x , h y is obtained througha double differentiation of ( . ). A first differentiation yields˙ h x = − l ˙ θ sin θ − l ( ˙ θ + ˙ θ ) sin ( θ + θ ) ˙ h y = l ˙ θ cos θ + l ( ˙ θ + ˙ θ ) cos ( θ + θ ) And then¨ h x = − h y ¨ θ − l sin ( θ + θ ) ¨ θ − φ x ( θ , ˙ θ ) ( . a)¨ h y = h x ¨ θ + l cos ( θ + θ ) ¨ θ − φ y ( θ , ˙ θ ) ( . b)with φ x ( θ , ˙ θ ) = l ˙ θ cos θ + l ( ˙ θ + ˙ θ ) cos ( θ + θ ) φ y ( θ , ˙ θ ) = l ˙ θ sin θ + l ( ˙ θ + ˙ θ ) sin ( θ + θ ) Thus, one gets (cid:32) ¨ h x ¨ h y (cid:33) = (cid:32) − h y − l sin ( θ + θ ) h x l cos ( θ + θ ) (cid:33) (cid:32) ¨ θ ¨ θ (cid:33) − (cid:32) φ x φ y (cid:33) Or, in other terms¨ h = H ¨ θ − φ . . flatness appls . for simplistic motor control 27 With the following notations H = (cid:32) − h y − l sin ( θ + θ ) h x l cos ( θ + θ ) (cid:33) , h = (cid:32) h x h y (cid:33) ( . )And, using ( . ), one gets the dynamics in the end effector h :¨ h = − H M − (cid:0) C + G − T (cid:1) − φ ( . ) . . Trajectory tracking of the two link arm
The system ( . ) is differentially flat, with flat output h x , h y . Indeedequations ( . ) yield the expressions of θ and θ in terms of h x , h y ,and T is given by: T = C + G + MH − (cid:0) ¨ h + φ (cid:1) Thus, considering a reference trajectory h xr ( t ) , h yr ( t ) , the so-calledopen loop control T r is given by: T r = C r + G r + MH − r (cid:0) ¨ h r + φ r (cid:1) ( . )with C r = (cid:32) − m l θ r r sin θ r − m l ˙ θ r ˙ θ r r sin θ r m l ˙ θ r r sin θ r (cid:33) G r = (cid:32) ( m l + m r ) g sin θ r + m r g sin ( θ r + θ r ) m r g sin ( θ r + θ r ) (cid:33) H r = (cid:32) − h yr − l sin ( θ r + θ r ) h xr l cos ( θ r + θ r ) (cid:33) φ r = (cid:32) l ˙ θ r cos θ r + l ( ˙ θ r + ˙ θ r ) cos ( θ r + θ r ) l ˙ θ r sin θ r + l ( ˙ θ r + ˙ θ r ) sin ( θ r + θ r ) (cid:33) θ r = arctan (cid:18) h yr h xr (cid:19) − arctan (cid:18) l sin θ r l + l cos θ r (cid:19) θ r = arctan (cid:32) ± (cid:112) − ¯ h r ¯ h r (cid:33) ¯ h r = h xr + h yr − l − l l l . flatness applications & extensions Then, the feedback control law ensuring tracking of the referencetrajectory h xr ( t ) , h yr ( t ) is given by: T = C + G + MH − (cid:16) φ + ¨ h r − Λ h e h − Λ h ˙ e h (cid:17) ( . )with Λ h = (cid:32) λ h λ h (cid:33) , Λ h = (cid:32) λ h λ h (cid:33) where the λ ijh are suitably chosen reals such that the closed loop errorequation in e h is exponentially stable (it is thus sifficient to choosethese as strictly positive reals).Note that the difference between the previous tracking control lawand the feedforward one given in ( . ) is of the form: T − T r = C − C r + G − G r + MH − (cid:0) φ + ¨ h r (cid:1) − MH − r (cid:0) φ r + ¨ h r (cid:1) − MH − (cid:16) λ Th e h + µ Th ˙ e h (cid:17) ( . )which tends to zero when e h itself tends to zero.In ( . ), one needs to compute H − (the matrix H being definedin ( . )), which requires the determinant ∆ H ∆ H = l (cid:0) h x sin ( θ + θ ) − h y cos ( θ + θ ) (cid:1) to be non zero. From ( . ), we get ∆ H = l (cid:0) h x sin θ − h y cos θ (cid:1) Thus, when ∆ H =
0, we gettan ( θ + θ ) = tan θ Or, what is the same θ =
0, or θ = π The first case yields the following end effector coordinates: h x = ( l + l ) cos θ h y = ( l + l ) sin θ . . flatness appls . for simplistic motor control 29 Thus, the end effector with coordinates h x , h y remains on a circlecentered at the origin and with radius l + l , which corresponds tothe arm being fully extended. The second case ( θ = π ) yields theend effector coordinates: h x = ( l − l ) cos θ h y = ( l − l ) sin θ and the end effector with coordinates h x , h y remains on a circle cen-tered at the origin and with radius l − l , wich corresponds to thearm fully folded.When designing a reference trajectory, we shall avoid these twocases. Let us consider h yr ( t ) = h y f − h yi (cid:104) + tanh (cid:16) γ ( h xr ( t ) − h x ) (cid:17)(cid:105) ( . a) h xr ( t ) = ( h x f − h xi ) tT + h xi ( . b)for t ∈ [ T ] , and for example: h xi = ( l + l ) , h x f = . a) h yi = l + l , h y f = − l ( . b)The trajectory tracking is illustrated in Figure . a, and the associ-ated animation in Figure . b. The corresponding control laws are x (m) h y ( m ) (a) End effector tracking: h y versus h x . x (m) h y ( m ) (b) Animated tracking of the two link arm. Figure . : Trajectory tracking of a two link robot arm. In red and dashedthe reference trajectory; in blue and solid, the actual (simulated)trajectory. shown in Figures . a and . b. The tracking errors are plotted in . flatness applications & extensions T , T r ( N ) (a) Tracking feedback control T . T , T r ( N ) (b) Tracking feedback control T . Figure . : Two link arm tracking feedback control laws. In red and dashedthe open loop (feedforward) law; in blue and solid, the feedbacklaw. −3 time (s) h x − h x r ( m ) (a) Tracking feedback error h x − h xr . h y − h y r ( m ) (b) Tracking feedback error h y − h yr . Figure . : Two link arm tracking feedback errors. Figures . a and . b.The previous results can be applied to more complete and less sim-plistic models, as in e.g. (Friston, ); see also (Richardson et al., ). . Several extensions of differentially flat systems can be envisioned.The recent paper (Aschenbrenner et al., ) reviews some of themost interesting ones. A Liouvillian closed structure D will con-tain all the solutions of first order linear differential equations, an . extended Liouvillian one will contain all the solutions of linear dif-ferential equations and an existentially closed one all the solutions ofalgebraic differential equations. We shall give some rather elementarydefinitions below (see Appendix C. , p. ). Definition The system ˙ x = f ( x , u ) with x ( t ) ∈ R n and u ( t ) ∈ R m is called Liouvillian (resp. extendedLiouvillian, existentially closed ) if there exists a set of variables, calleda
Liouvillian (resp. extended Liouvillian, existentially closed ) output y = ( y , . . . , y m ) solution ofH ( y , ˙ y , . . . , y ( r y ) , x , u , ˙ u , . . . , u ( r u ) ) = r y , r u ∈ N ( . ) with H linear of first order (resp. linear, polynomial) in its variables, suchthat x = A ( y , ˙ y , . . . , y ( ρ x ) ) u = B ( y , ˙ y , . . . , y ( ρ u ) ) with q an integer, and such that the system equationsdAdt ( y , ˙ y , . . . , y ( q + ) ) = f ( A ( y , ˙ y , . . . , y ( q ) ) , B ( y , ˙ y , . . . , y ( q + ) )) are identically satisfied. . . . Jansen and Rit modelBrief recall of the model
Consider the Jansen and Rit model, as depicted in (Pinotsis et al., ):¨ ν + κ e ˙ ν + κ e ν = κ e m e (cid:0) w F ( ν ) + u (cid:1) ( . a)¨ ν + κ i ˙ ν + κ i ν = κ i m i w F ( ν ) ( . b)¨ ν + κ e ˙ ν + κ e ν = κ e m e (cid:0) w F ( ν ) + w F ( ν ) (cid:1) ( . c) y = ν ( . d) . flatness applications & extensions These equations respectively depict the following populations: excita-tory stellate, inhibitory and excitatory. The signification of the variousvariables and parameters are the following: ν i Expected depolarization in the i -th population w ij F ( ν j ) presynaptic input to the i -th population from the j th one F ( ν j ) Sigmoid function of the postsynaptic depolarization w ij Instrinsic connection strength between the populations j and jm i , m e Maximum postsynaptic responses κ e , κ i Rate constants of postsynaptic filtering u Exogenous input y Endogenous outputThe choice made in (Pinotsis et al., ) for the sigmoid function F is the logistic function F ( ν ) = + e − β ( ν − ν T ) whose derivative and inverse are: F (cid:48) = β F ( F − ) , and F − ( η ) = φ ( η ) = ν T + β ln ηη − d i , d e be the differential operators d i = d dt + κ i ddt + κ i = (cid:18) ddt + κ i (cid:19) d e = (cid:18) ddt + κ e (cid:19) Then, the previous model ( . ) can be written as d e ν = κ e m e (cid:0) w F ( ν ) + u (cid:1) ( . a) d i ν = κ i m i w F ( ν ) ( . b) d e ν = κ e m e (cid:0) w F ( ν ) + w F ( ν ) (cid:1) ( . c) y = ν ( . d) Differential flatness of the model
A flat output of the model ( . ) is ν . Indeed, after equation ( . b),one gets ν : ν = φ (cid:18) κ i m i w ( ¨ ν + κ i ˙ ν + κ i ν ) (cid:19) ( . ) . Then, after ( . c) w F ( ν ) = w F ( ν ) + κ e m e ( ¨ ν + κ e ˙ ν + κ e ν ) Hence the expression for ν : ν = φ (cid:20) w w F ( ν ) + κ e m e w ( ¨ ν + κ e ˙ ν + κ e ν ) (cid:21) ( . )And, using ( . a), the expression for u : u = − w F ( ν ) + κ e m e ( ¨ ν + κ e ˙ ν + κ e ν ) ( . )The Figure . below outlines the compartmental like model under-lying the model ( . ). The bold arrows in this Figure enables one, u ν ν ν y input w F ( . ) L e ( . ) w F ( . ) w F ( . ) L e ( . ) w F ( . ) L i ( . ) ouput Figure . : The Jansen and Rit Model. by reversing the arrows, to reveal the differential flatness character ofthe model: ν is obtained from ν by reversing the arrow ( ν ) → ( ν ) (yielding equation ( . )) ; then, ν is obtained from ν (and ν ) by re-versing the arrow ( ν ) → ( ν ) (yielding equation ( . )) ; finally u isobtained from ν (and ν ) by reversing the arrow ( u ) → ( ν ) (yieldingequation ( . )). . flatness applications & extensions Extended Liouvillian character of the model
Since the output of interest considered in (Pinotsis et al., ) is ν , we can investigate how the model can be parametrized by thisvariable. The variable ν can be obtained from ν by integrating thedifferential equation ( . b) in ν (which is linear in this variable).Indeed, ( . b) can be rewritten as: ddt (cid:32) ν ˙ ν (cid:33) = (cid:32) − κ i − κ i (cid:33) (cid:32) ν ˙ ν (cid:33) + (cid:32) κ i m i w F ( ν ) (cid:33) or, in matrix form˙ V = AV + U , with V = (cid:32) ν ˙ ν (cid:33) , A = (cid:32) − κ i − κ i (cid:33) , U = (cid:32) κ i m i w F ( ν ) (cid:33) The general solution of this last equation is well known to be V = V ( ) e At + (cid:90) t e A ( t − τ ) U ( τ ) d τ One has e At = (cid:32) ( + κ i t ) e − κ i t te − κ i t − κ i e − κ i t ( − κ i t ) e − κ i t (cid:33) Thus, ν is given by ν = ν (cid:0) ( + κ i t ) e − κ i t (cid:1) + ˙ ν te − κ i t + κ i m i w (cid:90) t ( t − τ ) e − κ i ( t − τ ) F ( ν ( τ )) d τ ( . )where ν = ν ( ) , ˙ ν = ˙ ν ( ) . Then, the two other variables areobtained as in ( . )–( . ): ν = φ (cid:20) w w F ( ν ) + κ e m e w ( ¨ ν + κ e ˙ ν + κ e ν ) (cid:21) ( . a) u = − w F ( ν ) + κ e m e ( ¨ ν + κ e ˙ ν + κ e ν ) ( . b)Recalling d i , d e the differential operators d i = (cid:18) ddt + κ i (cid:19) d e = (cid:18) ddt + κ e (cid:19) . The previous equations ( . )–( . ) can be rewritten as: ν = d − i (cid:0) κ i m i w F ( ν ) (cid:1) ( . a) ν = φ (cid:20) w w F ( ν ) + κ e m e w d e ν (cid:21) ( . b) u = − w F ( ν ) + κ e m e d e ν ( . c)Thus, the model is extended Liouvillian and an extended Liouvillianoutput is ν . art IIN E U R A L F I E L D P O P U L A T I O N M O D E L S . . . General case model
We consider spatially distributed network models, such as the onesconsidered in Chapter of (Ermentrout and Terman, ), Subsec-tion . , p. , and Chapter , Subsection . . , p. or in Chapter of (Bressloff, ), Subsection . , p. (as well as Subsection . ,p. of (Bressloff, )).We can consider the so-called activity-based neural field model: τ sy ∂ν ( t , x ) ∂ t = − ν ( t , x ) + F (cid:0) I r w ( x ) ∗ x ν ( t , x ) (cid:1) h ( ν ( t , x )) + u ( t , x ) or, in a slightly more compact way τ sy ∂ν∂ t = − ν + F (cid:0) I r w ( x ) ∗ x ν (cid:1) h ( ν ) + u ( . )And a slightly more general case τ sy ∂ν∂ t = − ν + F (cid:18) I r (cid:90) Ω x w ( x − ξ ) ν ( t , ξ ) d ξ (cid:19) h ( ν ) + u ( . )Alternately, we can consider the slightly different model τ sy ∂ν∂ t = − ν + I r w ( x ) ∗ x F (cid:0) ν ( t , x ) (cid:1) h ( ν ) + u ( . )And its slight generalization τ sy ∂ν∂ t = − ν + I r (cid:90) Ω x w ( x − ξ ) F (cid:0) ν ( t , ξ ) (cid:1) d ξ h ( ν ) + u ( . ) . . Parameters and variables assumptions
The various parameters and functions staisfy the following: ◦ The spatial variable is three dimensional, i.e. Ω x ⊂ R . ◦ The parameter τ sy is a constant. ◦ If the synapses are saturating, then h ( s ) = − s , otherwise, h ( s ) = ◦ The function F is a sigmoid type function (see A. , p. ). ◦ The control u ( t , r ) has a spatial compact support Ω u . ◦ The neuron interaction strength function w ( r ) is symmetric, noneg-ative, integrates to over the whole line and is rapidly decayingat infinity: ∃ M ∈ R , ∃ α ∈ R + , ∀ x ∈ R s.t. (cid:107) x (cid:107) > M , (cid:107) w ( x ) (cid:107) < (cid:107) x − α (cid:107) . flatness applications & extensions Acronym Name Function w (Wdorg) Dirac at the origin δ (Wdnor) Dirac not at the origin δ x (Wsofd) Sum of Diracs ∑ Ni = a i δ x i (Wsexp) Single exponential e − ax H ( t ) (Wmexp) Multiple exponential ∑ Ni = e − a i x H ( t ) (Wgaus) Gaussian e − x / σ (Waexp) Absolute exponential e −| x | /2 (Wdosc) Decaying oscillatory e − b | x | ( b sin | x | + cos x ) (Wfhat) Flat hat shaped rect ( x / χ ) (Wmhat) Mexican hat Γ e − γ x − Γ e − γ x or e − c x / σ / σ − e − c x / σ / σ (Wwhat) Wizard hat ( )( − | x | ) e −| x | (Wcomp) Compact support w has compact support
Table . : Neuron interaction strength functions examples. . . Neuron interaction strength examples
Some typical examples of such w functions are shown in Table . .The graphical representations of some of the above mentionnedneuron interaction strength functions are given in Figures . and . . . . Simplification hypotheses
We shall make the following simplifying assumptions:(H ) We consider a spherical symmetric case, which boils down tothe unidimensional case where Ω x = [ a , b ] in r and r = (cid:107) x (cid:107) .(H ) We consider F equal to a heaviside: ∀ ξ ∈ R − , F ( ξ ) = ∀ ξ ∈ R + , F ( ξ ) = . −10 −5 0 5 1000.20.40.60.81 x w ( x ) (a) Absolute exponential. −10 −5 0 5 10−0.500.51 x w ( x ) (b) Decaying oscillatory with b = Figure . : Examples of neuron interaction strength functions. −10 −5 0 5 10−0.200.20.40.60.8 x w ( x ) (a) Mexican hat through gaussian differ-ence with c = c = m = m = −10 −5 0 5 10−0.100.10.20.3 x w ( x ) (b) Wizard hat. Figure . : Examples of neuron interaction strength functions. (H ) We consider that the synapses are not saturating, hence h ( s ) = . . Pointwise neuron interaction strength modelsDirac at the origin case
We consider that w is a single spatial Dirac: w ( r ) = δ . flatness applications & extensions We then obtain the following system˙ ν ( t , r ) = − a ν ( t , r ) + δ ∗ r (cid:0) [ a , b ] ν ( t , r ) (cid:1) + u ( t , r ) Or, in other words˙ ν ( t , r ) = − a ν ( t , r ) + [ a , b ] ( r ) ν ( t , r ) + u ( t , r ) ( . ) Dirac not at the origin case
We consider that w is a single spatial Dirac: w ( t , r ) = δ r with r ∈ [ a , b ] . We then obtain the following system˙ ν ( t , r ) = − a ν ( t , r ) + δ r ∗ r (cid:0) [ a , b ] ν ( t , r ) (cid:1) + u ( t , r ) ( . )Or, in other words˙ ν ( t , r ) = − a ν ( t , r ) + [ a , b ] ( r − r ) ν ( t , r − r ) + u ( t , r ) ( . ) . . Exponential type neuron interaction strengthA single exponential
Consider that the neuron interaction strength w satisfies a linear dif-ferential equation: w (cid:48) ( r ) = − aw ( r ) The model is the following τ sy ∂ν ( t , r ) ∂ t = − ν ( t , r ) + I r w ( r ) ∗ r ν ( t , r ) + u ( t , r ) ( . )Hence, the convolution part is w ( r ) ∗ r ν ( t , r ) = I r (cid:18) τ sy ∂ν ( t , r ) ∂ t + ν ( t , r ) − u ( t , r ) (cid:19) And is spatial derivative is w (cid:48) ( r ) ∗ r ν ( t , r ) = I r (cid:18) τ sy ∂ ν ( t , r ) ∂ t ∂ r + ∂ν ( t , r ) ∂ r + ∂ u ( t , r ) ∂ r (cid:19) = − aw ( r ) ∗ r ν ( t , r )= − aI r (cid:18) τ sy ∂ν ( t , r ) ∂ t + ν ( t , r ) − u ( t , r ) (cid:19) Hence, ν satisfies the following partial differential equation: τ sy ∂ t ∂ r ν ( t , r ) = − ( a τ sy ∂ t + ∂ r + a ) ν ( t , r ) − ( ∂ r + a ) u ( t , r ) . A more general case
In (Coombes et al., ), Chapter , “PDE Methods for Two-DimensionalNeural Fields” by Carlo R. Laing, the case of neuron interactionstrength w with a rational Fourier transform is considered: F ( w )( ξ ) = (cid:101) w ( ξ ) = (cid:90) + ∞ − ∞ w ( τ ) e − j ξτ d τ = p ( ξ ) q ( ξ ) with p and q two polynomials. Then, equation ( . ) with h ( s ) = τ sy ∂ν ( t , x ) ∂ t = − ν ( t , x ) + I r w ( x ) ∗ x F (cid:0) ν ( t , x ) (cid:1) + u ( t , x ) becomes ( τ sy ∂ t + ) ˜ ν ( t , ξ ) = I r (cid:101) w ( ξ ) (cid:101) F ( ν )( t , ξ ) + ˜ u ( t , ξ )= I r p ( ξ ) q ( ξ ) (cid:101) F ( ν )( t , ξ ) + ˜ u ( t , ξ ) Thus, through multiplication by q ( ξ ) : ( τ sy ∂ t + ) q ( ξ ) ˜ ν ( t , ξ ) = I r p ( ξ ) (cid:101) F ( ν )( t , ξ )+ q ( ξ ) ˜ u ( t , ξ ) which yields, in the spatial domain: ( τ sy ∂ t + ) q ( ∂ x ) ν ( t , x ) = I r p ( ∂ x ) F ( ν )( t , x )+ q ( ∂ x ) u ( t , x ) The Fourier tranforms of some of the neuron interaction strengthfunctions of Table . , p. are shown in Table . . . Let us consider, after (Pinotsis et al., ), the following neural fieldJansen and Rit model ¨ ν + κ e ˙ ν + κ e ν = κ e m e (cid:0) µ + u (cid:1) ( . a)¨ ν + κ i ˙ ν + κ i ν = κ i m i µ ( . b)¨ ν + κ e ˙ ν + κ e ν = κ e m e µ ( . c) ∂ t µ − σ ∂ x µ + σβ ∂ t µ + σ β µ = φ ( ν ) ( . d) ∂ t µ − σ ∂ x µ + σβ ∂ t µ + σ β µ = φ ( ν ) ( . e) ∂ t µ − σ ∂ x µ + σβ ∂ t µ + σ β µ = ψ ( ν , ν ) ( . f) . flatness applications & extensions Name Function/Distribution Fourier transform
Dirac at theorigin δ Dirac not atthe origin δ x e − j ξ x Sum ofDiracs N ∑ i = a i δ x i N ∑ i = a i e − j ξ x i Flat hatshaped rect (cid:18) x χ (cid:19) sinc Gaussian e − x σ σ √ e − σ ξ Absoluteexponential e − a | x | aa + ξ Decayingoscillatory e − b | x | ( b sin | x | + cos x ) b ( b + ) ξ + ( b − ) ξ + ( b + ) Mexican hat Γ e − γ | x | − Γ e − γ | x | Γ γ ( γ + ξ ) − Γ γ ( γ + ξ )( γ + ξ )( γ + ξ ) Wizard hat ( − | x | ) e −| x | ξ ( + ξ ) Table . : Fourier tranforms of some neuron interaction strength functions. . with φ ij ( ν j ) = α ij (cid:0) σ F ( ν j ) + σ F (cid:48) ( ν j ) (cid:1) ψ ( ν , ν ) = α (cid:16) σ β (cid:0) F ( ν ) − F ( ν ) (cid:1) + σ (cid:0) F (cid:48) ( ν ) − F (cid:48) ( ν ) (cid:1)(cid:17) Let us recall the differential operators d i and d e and introduce D , D , D : d i = (cid:18) ddt + κ i (cid:19) d e = (cid:18) ddt + κ e (cid:19) D = ∂ t − σ ∂ x + σβ ∂ t + σ β D = ∂ t − σ ∂ x + σβ ∂ t + σ β , D = ∂ t − σ ∂ x + σβ ∂ t + σ β The model ( . ) can then be rewritten as d e ν = κ e m e (cid:0) µ + u (cid:1) ( . a) d i ν = κ i m i µ ( . b) d e ν = κ e m e µ ( . c) D µ = φ ( ν ) ( . d) D µ = φ ( ν ) ( . e) D µ = ψ ( ν , ν ) ( . f)The Figure . below outlines the compartmental like model underly-ing the model ( . ).The variable ν is a Liouvillian output. Indeed, µ is obtainedthrough ν using equation ( . b): µ = κ i m i d i ν ( . )Then ν is obtained via µ with the help of ( . e) ν = ψ ( D µ ) ( . )where ψ denotes the inverse function of α ( σ β F + σ F (cid:48) ) . After, µ is derived from ν through ( . c) µ = κ e m e d e ν ( . ) . flatness applications & extensions u ν ν ν µ µ µ input Legend(1) α (cid:0) σ β F ( . ) + σ F (cid:48) () (cid:1) (2) α (cid:0) σ β F ( . ) + σ F (cid:48) ( . ) (cid:1) (3) α (cid:0) σ β F ( . ) + σ F (cid:48) ( . ) (cid:1) (4) − α (cid:0) σ β F ( . ) + σ F (cid:48) ( . ) (cid:1) κ e m e κ i m i κ e m e d e d i d e (2) (3)(1) (4) D D D Figure . : A Jansen and Rit neural field model. The variable µ and ν yields ν through ( . f) ν = D µ + α (cid:16) σ β F ( ν ) + σ F (cid:48) ( ν ) (cid:17) ( . )where ψ is the inverse function of α ( σ β F + σ F (cid:48) ) . By integratingthe wave equation in the equation ( . d) we get µ from ν µ = α D − (cid:0) σ β F ( ν ) + σ F (cid:48) ( ν ) (cid:1) ( . )At last, the control input u can de derived through ( . a) from µ and ν u = κ e m e d e ν − µ ( . ) I B L I O G R A P H YM. Aschenbrenner, L. van den Dries, and J. van der Hoeven. Towardsa model theory for transseries. Notre Dame Journal of Formal Logic , : – , .M. Aschenbrenner and L. van den Dries. Liouville closed h-fields. Journal of Pure and Applied Algebra , : – , a.M. Aschenbrenner and L. van den Dries. Analyzable Functions andApplications , volume of Contemp. Math. , chapter Asymptotic dif-ferential algebra, pages – . Amer. Math. Soc., Providence, RI, b.P.C. Bressloff. Spatiotemporal dynamics of continuum neural fields. J. Phys. A: Math. Theor. , , .P.C. Bressloff. Waves in Neural Media – From Single Neurons to NeuralFields . Springer, New York, .R. Brette and W. Gerstner. Adaptive exponential integrate-and-firemodel as an effective description of neuronal activity.
J. Neurophys-iol. , : – , .S. Coombes, P. beim Graben, R. Potthast, and J. Wright. Neural Fields– Theory and Applications . Springer, Berlin, .P. Dayan and L. F. Abbott.
Theoretical Neuroscience: Computational andMathematical Modeling of Neural Systems . The MIT Press, Cambridge,Massachussets, .B. Dumitrescu.
Positive Trigonometric Polynomials and Signal ProcessingApplications . Springer, Dordrecht, The Netherlands, .G.B. Ermentrout and D.H. Terman.
Mathematical Foundations of Neu-roscience . Springer, New York, .M. Fliess, J. Lévine, P. Martin, and P. Rouchon. Flatness and defect ofnon-linear systems: introductory theory and applications.
Internat.J. Control , : – , . Bibliography
M. Fliess, H. Mounier, P. Rouchon, and J. Rudolph. Tracking controlof a vibrating string with an interior mass viewed as delay system.
ESAIM Control Optim. Calc. Var. , : – , .K. Friston. The free-energy principle: a unified brain theory? Nat.Rev. Neurosci. , : – , .K. Friston. A free energy principle for biological systems. Entropy , : – , .H. Haken. Brain Dynamics – An Introduction to Models and Simulation .Springer-Verlag, Berlin, .E.M. Izhikevich. Simple model of spiking neurons.
IEEE Transactionson Neural Networks , : – , .E.M. Izhikevich. Which model to use for cortical spiking neurons? IEEE Transactions on Neural Networks , : – , .E.M. Izhikevich. Hybrid spiking models. Phil. Trans. R. Soc. A , : – , .P. Martin and P. Rouchon. Systèmes plats : planification et suivi detrajectoires. In Actes des Journées Nationales de Calcul Formel , pages – , . URL http://jncf2008.loria.fr/jncf2008.pdf .K.D. Miller and F. Fumarola. Mathematical equivalence of two com-mon forms of firing-rate models of neural networks. Neural Comput , : – , .K. Friston and D.A. Pinotsis. Xx.D.A. Pinotsis, R.J. Moran, and K.J. Friston. Dynamic causal modelingwith neural fields. Neuroimage , : – , .M.J. Richardson, M.A. Riley, and K. Shockley, editors. Progress in Mo-tor Control – Neural, Computational and Dynamic Approaches , volume of Advances in experimental medicine and biology . Springer, Dor-drecht, .P. Rouchon. Motion planning, equivalence, infinite dimensional sys-tems.
Int. J. of Applied Mathematics and Computer Science , , .J. Touboul. Importance of the cutoff value in the quadratic adaptiveintegrate-and-fire model. Neural Comput. , : – , . ibliography J. van der Hoeven.
Transseries and real differential algebra , volume of Lecture Notes in Mathematics . Springer-Verlag, .F. Woittennek and H. Mounier. Controllability of networks of spa-tially one-dimensional second order PDE – an algebraic approach.
Siam J. Contr. , : – , . art IIIA P P E N D I X N O TAT I O N S , T R A N S F O R M S A N D S I G M O I DF U N C T I O N S a . ◦ The distribution H ( η ) is the Heaviside distribution: H ( η ) = η (cid:54)
01 if η > ◦ The function sinc ( η ) is the cardinal sine :sinc ( η ) = sin ( η ) η ◦ The distribution rect ( η ) is the rectangular pulse of width 1:rect ( η ) = | η | (cid:62)
121 if | η | < ◦ The boxcar distribution (rectangular pulse of width ρ and cen-tered on η ):rect (cid:18) t − t ρ (cid:19) = H ( η − η + ρ ) − H ( t − t − ρ ) ◦ The linear rectifier function is | η | + = v (cid:54) v when v >
534 notations , transforms and sigmoid functions a . Consider a function f ( t , x ) from R × Ω x to R , where Ω x ⊆ R . ◦ The function ˜ f ( t , ξ ) will designate the spatial Fourier transform of f , i.e.˜ f ( t , ξ ) = F ( f )( t , ξ ) = (cid:90) + ∞ − ∞ f ( t , x ) e − j ξ x dx ◦ The function ˆ f ( s , x ) will designate the temporal Laplace transform of f , i.e.ˆ f ( s , x ) = L ( f )( s , x ) = (cid:90) + ∞ − ∞ f ( t , x ) e − st dt a . The following functions F will be in particular used for the firingrate. Thus F ( ξ ) designates a spike rate and ξ a sitmulus intensity. Asigmoid function F : R → R is such that F ( ) = F (cid:48) ( ) = F ( ) = F (cid:48) ( ) = ∀ ξ ∈ R , ξ < F ( ξ ) = ∀ ξ ∈ R , ξ > F ( ξ ) = ), Section . , p. ; (Bressloff, ),pp. , , , , (Haken, ), pp. , ). ◦ The
Heaviside . F ( x ) = F H ( ξ − ξ ) ◦ The
Piecewise linear function. F ( x ) = ξ < ξ β ( ξ − ξ ) if ξ (cid:54) ξ < ξ + β ξ > ξ + β ◦ The
Logistic function. F ( x ) = + e − β ( x − x T ) . one has F (cid:48) = β F ( F − ) F − ( η ) = φ ( η ) = x T + β ln ηη − ◦ The
Traub Model . F ( ξ ) = + e − ξ − βα ◦ The
Hyperbolic tangent function. F ( ξ ) = F (cid:0) + tanh ( αξ ) (cid:1) ◦ The
Square root function. F ( ξ ) = F (cid:112) ξ − ξ T ◦ The
Noisy firing rate function. F ( ξ ) = (cid:115) ξ − ξ T − e − ( ξ − ξ T ) β Here, β is a measure of the noise, and when β tends to zero, thefunction approaches a pure square root model. ◦ The
Mean firing rate with flexible shape function (see (Coombeset al., ), p. ). F ( ξ ) = F m − F m (cid:0) + e √ ξ − µσ (cid:1) κ ◦ The
Naka-Rushton functions (alternately called Hill functions). F ( ξ ) = r ξ n ξ n + θ n if ξ (cid:54)
00 if ξ > r is the maximum spike rate and θ is the value of thestimulus intensity for which F reaches half its maximum. Theexponent n is a measure of the steepness of the F ( ξ ) curve. Typ-ical values matching experimental data range from 1.4 to 3.4.The function F ( ξ ) = − ξ n ξ n + θ = θ n ξ n + θ n is also used. , transforms and sigmoid functions ◦ The
Algebraic sigmoid function . F ( ξ ) = ξ (cid:112) + ξ This function has the inverse F ( η ) = η (cid:112) − η S O M E F L AT N E S S S I M P L E C R I T E R I AThere does not exist, at the time of this writing, a general criterion forchecking flatness, neither for building flat outputs in a constructivemanner. Nevertheless, some peculiar cases are to be noticed. b . Proposition Any static state feedback linearizable system is flat.
See below (Subsection B. , p. ) for a static state feedback lineariz-ability criterion for affine input systems. Proposition (Charlet, Levine and Marino, ) For systems with asingle input, dynamic feedback linearization implies static feedback lineariza-tion.
Proposition (Charlet, Levine and Marino, ) A dynamics affine inthe input with n states and n − inputs is flat as soon as it is controllable(strongly accessible). Recall that a dynamics is called affine in the input if it is of the form˙ x = f ( x ) + n − ∑ i = g i ( x ) u i A dynamics with x ∈ X ⊆ R n is strongly accessible if, for all x ∈ X ,there exists a T > R ( x ) (cid:54) = ∅ where int S denotes the interior of the set S and R ( x ) is the reachableset of x .
578 some flatness simple criteria b . Proposition (Ruled variety criterion, Rouchon, ) Suppose the dy-namics ˙ x = f ( x , u ) is flat. The projection of the sub variety p = f ( x , u ) inthe ( p , u ) -space ( x is here a parameter) onto the p -space is a ruled varietyfor all x . This criterion means that the elimination of u from the n equations˙ x = f ( x , u ) yields n − m equations F ( x , ˙ x ) = ( x , p ) such that F ( x , p ) =
0, there exists a ∈ R n , a (cid:54) = ∀ λ ∈ R , F ( x , p + λ a ) = F ( x , p ) is thus ruled since it contains the line passingthrough p with direction a . b . Consider an affine input system˙ x = f ( x ) + n − ∑ i = g i ( x ) u i = f ( x ) + g ( x ) u where f , g i are smooth vector fields on a domain D ⊂ R n , x ∈ D , u ∈ R m . b . . Brief recall of differential geometry notions
Definition Let r (cid:62) be an integer. A C r vector field on R n is a map-ping f : D → R n of class C r from an open set D ⊂ R n to R n . A smoothvector field is a mapping f : D → R n of class C ∞ . Let h ( x ) be a smooth vector field on a domain D ∈ R n . The Liederivative of h along f , denoted as L f h ( x ) can be defined (in localcoordinates) as ∂ h ( x ) ∂ x f ( x ) = n ∑ i = ∂ h ( x ) ∂ x i f ( x ) since it is a smooth vector field, a Lie derivative operator can be ap-plied to it. Set L if = L f L i − f . The
Lie bracket of f and g can be defined (in local coordinates) as [ f , g ]( x ) = ∂ g ∂ x f ( x ) − ∂ f ∂ x g ( x ) Iterated lie brackets are denoted as ad if g :ad f g = [ f , g ] ad if g = [ f , ad i − f g ] Let f , . . . , f η be some vector fields on D ⊂ R n . The distribution ∆ spanned the vector fields f , . . . , f η is the collection of vector spaces ∆ ( x ) = span R n (cid:8) f ( x ) , f ( x ) , . . . , f η ( x ) (cid:9) for all x ∈ D . We denote ∆ = span R n (cid:8) f , f , . . . , f η (cid:9) A distribution ∆ is involutive if ∀ g , g ∈ ∆ , [ g , g ] ∈ ∆ Proposition The system with dyanmics ˙ x = f ( x ) + g ( x ) u is static statefeedback linearizable if, and only if, there is a domain D ⊂ D such that thefollowing two conditions are staisfied: . The matrix (cid:2) g , ad f g , . . . , ad n − f g (cid:3) has rank n for all x ∈ D . . The distribution (cid:8) g , ad f g , . . . , ad n − f g (cid:9) is involutive in D . P R E C I S E D E F I N I T I O N S F O R E X T E N S I O N S O FD I F F E R E N T I A L F L AT N E S S c . , dynamics and differential flatnessc . . Basic definitions from differential algebra
Definition An ordinary differential field k, is a field on which a map-ping d : k → k is defined, satisfying the natural properties with respect toaddition and product, i.e., for any x , z ∈ k, d ( x + z ) = d ( x ) + d ( z ) d ( xz ) = d ( x ) z + x d ( z ) Definition Let K be a field. A subfield of K is a subset k of K that isclosed under the field operations of K and under taking inverses in K. Inother words, k is a field with respect to the field operations inherited from K.The larger field K is then said to be an extension field of k, denoted as K / k. Definition Let k and K be differential fields with differential operators d k and d K respectively. Then, K is a differential extension field of k if K isan extension field of k and ∀ x ∈ k , d k ( x ) = d K ( x ) Let S be a subset of K. We shall denote by k (cid:104) S (cid:105) the differential subfield of Kgenerated by k and S. c . . Algebraic and transcendental extensions
All fields are assumed to be of characteristic zero. Assume also thatthe differential field extension K / k is finitely generated , i.e., there existsa finite subset S ⊂ K such that K = k (cid:104) S (cid:105) .
612 definitions for extensions of differential flatness
Definitions An element a of K is said to be differentially algebraic overk if it satisfies an algebraic differential equation with coefficients in k: thereexists a non-zero polynomial P over k, in several indeterminates, such thatP ( a , ˙ a , . . . , a ( ν ) ) = It is said to be differentially transcendental over k if it is not differentiallyalgebraic.The extension K / k is said to be differentially algebraic if any element ofK is differentially algebraic over k. An extension which is not differentiallyalgebraic is said to be differentially transcendental . c . . Nonlinear systems and flatness
Definition Let k be a given differential ground field. A (nonlinear) sys-tem is a finitely generated differential extension K / k. Definition A nonlinear system K / k is called differentially flat if thereexists a finite family y = ( y , . . . , y m ) of elements of an algebraic extensionL of K such that the extension L / k (cid:104) y (cid:105) is (non differentially) algebraic. Sucha family is called a flat output . c . - fields , liouvillian and existential closedness The paper (Aschenbrenner et al., ) reviews some of the most inter-esting notions for extending the differential flatness notion. One canalso see (Aschenbrenner and van den Dries, a,b, van der Hoeven, ) for related material.
Definition An H -field is an ordered differential field K whose naturaldominance relation (cid:52) satifies the following two conditions, for all z ∈ K:(H ) If z (cid:31) , then ˙ z / z > (H ) If z (cid:52) , then z − c ≺ , for some c ∈ C, the field of constants of K.
Remark In more usual terms, the dominance relations can be explicitedas follows, for real valued functions: ◦ f (cid:52) g ⇔ lim t → ∞ f ( t ) g ( t ) ∈ R ◦ f ≺ g ⇔ lim t → ∞ f ( t ) g ( t ) = . - fields , liouvillian and existential closedness 63 Definition An H-field K is
Liouville closed if it is real closed and anyequation ˙ z + az = b with a , b ∈ K has a non zero solution in K.
Definition An H-field K is existentially closed if every finite systemof algebraic differential equations over K in several unkowns with a solutionin an H-field extension of K has a solution in K.
T W O L I N K A R M I N V E R S E K I N E M AT I C SThe position of the end effector (the wrist) of a two link arm is givenby h x = l cos θ + l cos ( θ + θ ) (D. a) h y = l sin θ + l sin ( θ + θ ) (D. b)where h x , h y are the coordinates of the end effector. By summing thesquare of the two preceding equations, one obtains h x + h y = l + l + l l (cid:2) cos θ cos ( θ + θ ) + sin θ sin ( θ + θ ) (cid:3) = l + l + l l cos θ Then cos θ = h x + h y − l − l l l or θ = arctan (cid:18) sin θ cos θ (cid:19) = arctan (cid:32) ±√ − cos θ cos θ (cid:33) = arctan (cid:32) ± √ − ¯ h ¯ h (cid:33) ¯ h = h x + h y − l − l l l Setting k = l + l cos θ , k = l sin θ
656 two link arm inverse kinematics one has h x = k cos θ − k sin θ h y = k sin θ + k cos θ Then ρ = (cid:113) k + k , γ = arctan (cid:18) k k (cid:19) wherefrom h x = ρ cos γ cos θ − ρ sin γ sin θ h y = ρ cos γ sin θ + ρ sin γ cos θ or, what is the same h x ρ = cos ( γ + θ ) , h y ρ = sin ( γ + θ ) Then θ + γ = arctan (cid:18) h y h x (cid:19) and, finally θ = arctan (cid:18) h y h x (cid:19) − arctan (cid:18) l sin θ l + l cos θ (cid:19) T W O L I N K A R M C O D E L I S T I N G e . The listing beginning on the next page is a matlab code of the dif-ferential flatness based control of the two link arm example whosemodel is depicted in ( . ) and tracking feedback law in ( . ).
678 two link arm code listing % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Simulation of a two joint arm flatness based tracking% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function mainTwoJointArmFlatness ()% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Nomenclature% P: physical parameters% R: reference trajectories and conrol laws% U: complete control laws% G: gains % S: simulation scenario% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear all% %% PHYSICAL parameters% Fixed physical parametersP.g = 9.81; % gravity constantP. l1 = 0.3384; % (m) length of lower arm partP. l2 = 0.4554; % (m) length of upper arm part P. r1 = 0.1692; % (m) distance to center of massP. r2 = 0.2277; % (m) distance to center of massP. m1 = 2.10; % ( kg ) mass of lower arm partP. m2 = 1.65; % ( kg ) mass of upper arm partP. J1 = 0.025; % ( kg .m ^2) inertia of lower arm part P. J2 = 0.075; % ( kg .m ^2) inertia of upper arm part%% REFERENCE trajectory% A tanh oneR. hxi = 0.8*( P. l1 +P. l2 ); R. hxf = 0; R. hyi = P. l1 + 0.1* P. l2 ; R. hyf = - 0.1* P. l1 ;R. stiffness = 9; R. hxR = 0.5* R. hxi ;%% SIMULATION valuesS. tini = 0; S. tend = 10; % intial and ← (cid:45) final simulation times S. relTol = 1e -3; S. absTol = 1e -6; % relative and ← (cid:45) absolute simul tols% Intial state valuestvirt = [S. tini :0.01: S. tend ] ';strRef = inverseKinematics ( tvirt , P , R , S);S. theta10 = strRef . theta1r (1) + 0.05*( strRef . ← (cid:45) theta1r ( end ) - strRef . theta1r (1) ); S. dotTheta10 = strRef . dotTheta1r (1) + 0.05*( strRef . ← (cid:45) dotTheta1r ( end ) - strRef . dotTheta1r (1) );S. theta20 = strRef . theta2r (1) + 0.05*( strRef . ← (cid:45) theta1r ( end ) - strRef . theta1r (1) );S. dotTheta20 = strRef . dotTheta2r (1) + 0.05*( strRef . ← (cid:45) dotTheta2r ( end ) - strRef . dotTheta2r (1) ); . %% Default FEEDBACK GAINS % flatness based gains% for (s+ s1 )(s+ s2 ) = s ^2 + ( s1 + s2 )*s + s1 * s2sTh1 = 5; sDotTh1 = 2* sTh1 ;sTh2 = 6; sDotTh2 = 2* sTh2 ;K. kp1 = sTh1 * sDotTh1 ; K. kp2 = sTh2 * sDotTh2 ; K. kd1 = sTh1 + sDotTh1 ; K. kd2 = sTh2 + sDotTh2 ;K. Kp = [K. kp1 0; 0 K. kp2 ];K. Kd = [K. kd1 0; 0 K. kd2 ];% Saving option S. forSaving = ' yes ';% Simulateoptions = odeset ( ' RelTol ', S. relTol , ' AbsTol ', S. ← (cid:45) absTol );[ ts Xs ] = ode23tb ( @dynTwoJointArm , [S. tini S. tend ], ← (cid:45) ... [S. theta10 S. dotTheta10 S. theta20 S. ← (cid:45) dotTheta20 ] ' ,...options , P , R , K , S);S. ts = ts ; S. Xs = Xs ;plotVariables (P , R , K , S); end % of main ()function P = coriolisGravity (P , theta1 , dotTheta1 , theta2 ← (cid:45) , dotTheta2 )% Masses and intertia gathering m1 = P. m1 ; m2 = P. m2 ; l1 = P. l1 ; l2 = P. l2 ;r1 = P. r1 ; r2 = P. r2 ; J1 = P. J1 ; J2 = P. J2 ;% Coriolis and gravity terms computationsM11 = J1 + J2 + ( m1 * r1 ^2) + ( m2 *(( l1 ^2) + ( r2 ^2) + ← (cid:45) (2* l1 * r2 * cos ( theta2 ))));M12 = J2 + ( m2 *(( r2 ^2) + ( l1 * r2 * cos ( theta2 )))); M21 = M12 ;M22 = J2 + ( m2 * r2 ^2) ;M = [ M11 , M12 ; M21 , M22 ];C1 = -( m2 * l1 * r2 * dotTheta2 ^2* sin ( theta2 )) -...(2* m2 * l1 * r2 * dotTheta1 * dotTheta2 * sin ( theta2 )); C2 = m2 * l1 * dotTheta1 ^2* r2 * sin ( theta2 );C = [C1 , C2 ] ';G1 = (P.g* sin ( theta1 ) *(( m2 * l1 ) +( m1 * r1 ))) + (P.g* m2 * ← (cid:45) r2 * sin ( theta1 + theta2 ));G2 = P.g* m2 * r2 * sin ( theta1 + theta2 );G = [G1 , G2 ] '; P.C = C; P.G = G; P.M = M;end % %%%%%%%%%%% Dynamics % %%%%%%%%%%function [ dotX ] = dynTwoJointArm (t , X , P , R , K , S)persistent count ;if (t <= 0) count = 0; end ; if ( round (t) >= count )count = count + 1;disp ( sprintf ( ' time t = %f\n ', t));end ;% Variable gathering dotX = zeros (4 ,1) ;theta1 = X (1 ,:) ; dotTheta1 = X (2 ,:) ;theta2 = X (3 ,:) ; dotTheta2 = X (4 ,:) ;% Coriolis and gravity terms computationsP = coriolisGravity (P , theta1 , dotTheta1 , theta2 , ← (cid:45) dotTheta2 ); % Control law computationT = twoLinkFlatCtrlLaw (t , X , P , R , K , S);% Two link arm DYNAMICSM = P.M; C = P.C; G = P.G;ddotTheta = inv (M) * (-C - G + T); % return the derivative of the statedotX = [ dotTheta1 ddotTheta (1) dotTheta2 ddotTheta ← (cid:45) (2) ] ';end% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Flatness based tracking control law% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function [T] = twoLinkFlatCtrlLaw (t , X , P , R , K , S)% Variable gatheringdotX = zeros (4 ,1) ; theta1 = X (1 ,:) ; dotTheta1 = X (2 ,:) ;theta2 = X (3 ,:) ; dotTheta2 = X (4 ,:) ;M = P.M; C = P.C; G = P.G; l1 = P. l1 ; l2 = P. l2 ;% reference trajectories[ strRefHx strRefHy ] = tanhRefTraj (t , P , R , S); hxr = strRefHx .v; hyr = strRefHy .v;dotHxr = strRefHx . d1 ; dotHyr = strRefHy . d1 ;ddotHxr = strRefHx . d2 ; ddotHyr = strRefHy . d2 ;% Intermediary computations[ Hinv phi ] = computeHinvPhi ( theta1 , theta2 , ← (cid:45) dotTheta1 , dotTheta2 , P); hx = l1 * cos ( theta1 ) + l2 * cos ( theta1 + theta2 );hy = l1 * sin ( theta1 ) + l2 * sin ( theta1 + theta2 );dotHx = -l1 * dotTheta1 * sin ( theta1 ) - l2 *( dotTheta1 + ← (cid:45) dotTheta2 )* sin ( theta1 + theta2 ); . dotHy = l1 * dotTheta1 * cos ( theta1 ) + l2 *( dotTheta1 + ← (cid:45) dotTheta2 )* cos ( theta1 + theta2 );% errors computation ehx = hx - hxr ; ehy = hy - hyr ;dotEHx = dotHx - dotHxr ; dotEHy = dotHy - dotHyr ;% Vector computationvectDdotHr = [ ddotHxr ddotHyr ] ';vectEH = [ ehx ehy ] '; vectDotEH = [ dotEHx dotEHy ] ';% control law computationT = C + G + M* Hinv *( phi + vectDdotHr -...K. Kd * vectDotEH - K. Kp * ← (cid:45) vectEH );end % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Compute phi and the inverse of H% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function [ matHinv vectPhi ] = computeHinvPhi ( theta1 , ← (cid:45) theta2 , dotTheta1 , dotTheta2 , P) l1 = P. l1 ; l2 = P. l2 ;hx = l1 * cos ( theta1 ) + l2 * cos ( theta1 + theta2 );hy = l1 * sin ( theta1 ) + l2 * sin ( theta1 + theta2 );H = [- hy -l2 * sin ( theta1 + theta2 );hx l2 * cos ( theta1 + theta2 ) ]; % matHinv = [- l2 * cos ( theta1 + theta2 ) -l2 * sin ( theta1 + ← (cid:45) theta2 );% hx hy ← (cid:45) ]...% / ( l2 *( hy * cos ( theta1 + theta2 ) -hx * sin ( ← (cid:45) theta1 + theta2 )));matHinv = inv (H);vectPhi = [ l1 * dotTheta1 ^2* cos ( theta1 ) + l2 *( ← (cid:45) dotTheta1 + dotTheta2 ) ^2* cos ( theta1 + theta2 ); l1 * dotTheta1 ^2* sin ( theta1 ) + l2 *( ← (cid:45) dotTheta1 + dotTheta2 ) ^2* sin ( theta1 + ← (cid:45) theta2 ) ];end% %%%%%%%%%%%%%%%%%%%%%%% Reference trajectory % %%%%%%%%%%%%%%%%%%%%%%function [ strRefHx strRefHy ] = tanhRefTraj (t , P , R , S)hx = (( R. hxf - R. hxi ) .* t ) ./( S. tend ) + R. hxi ;strRefHx .v = hx ;strRefHx . d1 = (( R. hxf - R. hxi ) ./ S. tend ) .* ones ( length ( hx ← (cid:45) ) ,1) ; strRefHx . d2 = zeros ( length ( hx ) ,1) ;strRefTanh = tanhTr (hx , R. stiffness , R. hyi , R. hyf , R. ← (cid:45) hxR ); strRefHy .v = strRefTanh .v;strRefHy . d1 = strRefTanh . d1 .* strRefHx . d1 ;strRefHy . d2 = strRefTanh . d2 .* strRefHx . d1 .^2; end% %%%%%%%%%%%%%%%%%%%%% Inverse kinematics% %%%%%%%%%%%%%%%%%%%% function strRef = inverseKinematics (t , P , R , S)l1 = P. l1 ; l2 = P. l2 ;[ strRefHx strRefHy ] = tanhRefTraj (t , P , R , S);hxr = strRefHx .v; hyr = strRefHy .v;dothxr = strRefHx . d1 ; dothyr = strRefHy . d1 ; ddothxr = strRefHx . d2 ; ddothyr = strRefHy . d2 ;cosTheta2r = ( hxr .^2 + hyr .^2 - l1 .^2 - l2 .^2) ./ ← (cid:45) (2.* l1 .* l2 );sinTheta2r = sqrt (1 - cosTheta2r .^2) ;theta2r = unwrap ( atan2 ( sqrt (1 - cosTheta2r .^2) , ← (cid:45) cosTheta2r ));theta1r = unwrap ( atan2 ( hyr , hxr )) -... unwrap ( atan2 ( l2 .* sinTheta2r ,( l1 + l2 .* ← (cid:45) cosTheta2r )));dotTheta2r = -( hxr .* dothxr + hyr .* dothyr ) ./ ( l1 * l2 ← (cid:45) .* sinTheta2r );dotTheta1r = ( dothyr .* hxr - hyr .* dothxr ) ./ ( hxr .^2 ← (cid:45) + hyr .^2) -...( l2 .* dotTheta2r .*(1+ l1 .* cosTheta2r )) ← (cid:45) ./...( l1 ^2 + 2* l1 * l2 .* cosTheta2r + l2 ^2) ; ddotTheta2r = -( hxr .* ddothxr + dothxr .^2 + hyr .* ← (cid:45) ddothyr + dothyr .^2) ./...( l1 * l2 .* sinTheta2r ) + ( hxr .* dothxr + ← (cid:45) hyr .* dothyr ) .*...( cosTheta2r .* dotTheta2r ) ./ ( l1 * l2 .* ← (cid:45) sinTheta2r .^2) ;ddotTheta1r = ( ddothyr .* hxr - hyr .* ddothxr ) ./ ( hxr ← (cid:45) .^2 + hyr .^2) -...2.*( dothyr .* hxr - hyr .* dothxr ) .*( hxr .* ← (cid:45) dothxr + hyr .* dothyr ) -... l2 .*( ddotTheta2r .*(1+ l1 .* cosTheta2r ) - ← (cid:45) ( l1 .* dotTheta2r .^2.* sinTheta2r )) ← (cid:45) ./...( l1 ^2 + 2* l1 * l2 .* cosTheta2r + l2 ^2) + ← (cid:45)
2* l1 * l2 .*( l2 .* dotTheta2r .*(1+ l1 ← (cid:45) .*...cosTheta2r )) .*( sinTheta2r .* dotTheta2r ) ← (cid:45) ./ ( l1 ^2+2* l1 * l2 .* cosTheta2r + l2 ← (cid:45) ^2) .^2;strRef . theta1r = theta1r ; strRef . theta2r ← (cid:45) = theta2r ; . strRef . dotTheta1r = dotTheta1r ; strRef . dotTheta2r ← (cid:45) = dotTheta2r ; strRef . ddotTheta1r = ddotTheta1r ; strRef . ddotTheta2r ← (cid:45) = ddotTheta2r ;end% %%%%%%%% plots % %%%%%%%function plotVariables (P , R , K , S)% % For plotspolice = ' Helvetica '; size = 24; lineWidth = 2; % reference statets = S. ts ; tr = ts ; Xs = S. Xs ;% Reference curve and inverse kinematics[ strRefHx strRefHy ] = tanhRefTraj (tr , P , R , S); hxr = strRefHx .v; hyr = strRefHy .v;dotHxr = strRefHx . d1 ; dotHyr = strRefHy . d1 ;ddotHxr = strRefHx . d2 ; ddotHyr = strRefHy . d2 ;strRefTh = inverseKinematics (tr , P , R , S);theta1r = strRefTh . theta1r ; theta2r = ← (cid:45) strRefTh . theta2r ; dotTheta1r = strRefTh . dotTheta1r ; dotTheta2r = ← (cid:45) strRefTh . dotTheta2r ;% simulated statetheta1 = Xs (: ,1) ; dotTheta1 = Xs (: ,2) ;theta2 = Xs (: ,3) ; dotTheta2 = Xs (: ,4) ;% Control laws computation T1r = []; T2r = []; T1 = []; T2 = [];XsT = Xs ';for i = 1: length ( ts )% open loop controlPr = coriolisGravity (P , theta1r (i) , dotTheta1r (i) ← (cid:45) , theta2r (i) , dotTheta2r (i)); Mr = Pr .M; Cr = Pr .C; Gr = Pr .G;[ Hinvr phir ] = computeHinvPhi ( theta1r (i) , theta2r ← (cid:45) (i) , dotTheta1r (i) , dotTheta2r (i) , P);DdotHr = [ ddotHxr (i) ddotHyr (i) ] ';Tr = Cr + Gr - Mr * Hinvr *( DdotHr + phir );T1r = [ T1r ; Tr (1) ]; T2r = [ T2r ; Tr (2) ]; % simulated closed loop controlP = coriolisGravity (P , theta1 (i) , dotTheta1 (i) , ← (cid:45) theta2 (i) , dotTheta2 (i));T = twoLinkFlatCtrlLaw ( ts (i) , XsT (: , i) , P , R , K , ← (cid:45) S);T1 = [ T1 ; T (1) ]; T2 = [ T2 ; T (2) ];end ; % Simulated curve : forward kinematics l1 = P. l1 ; l2 = P. l2 ;hx = l1 .* cos ( theta1 ) + l2 .* cos ( theta1 + theta2 ) ← (cid:45) ;hy = l1 .* sin ( theta1 ) + l2 .* sin ( theta1 + theta2 ) ← (cid:45) ; figure (1) ;% hx hy plotsubplot (2 , 2, 1) ;plot (tr , hxr , 'r ', ts , hx , 'b ', ' LineWidth ', ← (cid:45) lineWidth ); grid ;xlabel ( ' time (s) ', ' FontName ', police , ' FontSize ', ← (cid:45) size ); ylabel ( 'h_x , h_ { xr } ( deg ) ', ' FontName ', police , ' ← (cid:45) FontSize ', size );title ( ' red h_ { xr } ; blue h_x ' ,...' FontName ', police , ' FontSize ', size , ' ← (cid:45) FontWeight ',' bold ');% theta2 plotsubplot (2 , 2, 2) ; plot (tr , hyr , 'r ', ts , hy , 'b ', ' LineWidth ', ← (cid:45) lineWidth ); grid ;xlabel ( ' time (s) ', ' FontName ', police , ' FontSize ', ← (cid:45) size );ylabel ( 'h_y , h_ { yr } ( deg ) ', ' FontName ', police , ' ← (cid:45) FontSize ', size );title ( ' red h_ { yr } ; blue h_y ' ,...' FontName ', police , ' FontSize ', size , ' ← (cid:45) FontWeight ',' bold '); subplot (2 , 2, 3) ;plot (tr , theta1r .*(180/ pi ) , 'r ', ts , theta1 .*(180/ pi ) ← (cid:45) , 'b ', ' LineWidth ', lineWidth ); grid ;xlabel ( ' time (s) ', ' FontName ', police , ' FontSize ', ← (cid:45) size );ylabel ( '\ theta_ {2} , \ theta_ {2 r} ( deg ) ', ' FontName ', ← (cid:45) police , ' FontSize ', size );title ( ' red \ theta_ {1 r} ; blue \ theta_1 ' ,... ' FontName ', police , ' FontSize ', size , ' ← (cid:45) FontWeight ',' bold ');subplot (2 , 2, 4) ;plot (tr , theta2r .*(180/ pi ) , 'r ', ts , theta2 .*(180/ pi ) ← (cid:45) , 'b ', ' LineWidth ', lineWidth ); grid ;xlabel ( ' time (s) ', ' FontName ', police , ' FontSize ', ← (cid:45) size );ylabel ( '\ theta_ {2} , \ theta_ {2 r} ( deg ) ', ' FontName ', ← (cid:45) police , ' FontSize ', size ); title ( ' red \ theta_ {2 r} ; blue \ theta_2 ' ,...' FontName ', police , ' FontSize ', size , ' ← (cid:45) FontWeight ',' bold '); . hFig2 = figure (2) ;% T1 and T1r plot set ( gca , ' FontName ', police , ' FontSize ', size );plot (tr , T1r , 'r -- ', ts , T1 , 'b - ', ' LineWidth ', ← (cid:45) lineWidth ); grid ;xlabel ( ' time (s) ', ' FontName ', police , ' FontSize ', ← (cid:45) size );ylabel ( 'T_1 , T_ {1 r} (N) ', ' FontName ', police , ' ← (cid:45) FontSize ', size );if ( strcmp (S. forSaving , ' yes ') ~= 1) title ( ' red T_ {1 r} ; blue T_1 ' ,...' FontName ', police , ' FontSize ', size , ' ← (cid:45) FontWeight ',' bold ');else print ( hFig2 , '- dpdf ', ' ../ GraphicsImages / ← (cid:45) twoLinkArmCtrl1 . pdf ');end hFig3 = figure (3) ;% T2 and T2r plotset ( gca , ' FontName ', police , ' FontSize ', size );plot (tr , T2r , 'r -- ', ts , T2 , 'b - ', ' LineWidth ', ← (cid:45) lineWidth ); grid ; xlabel ( ' time (s) ', ' FontName ', police , ' FontSize ', ← (cid:45) size );ylabel ( 'T_2 , T_ {2 r} (N) ', ' FontName ', police , ' ← (cid:45) FontSize ', size );if ( strcmp (S. forSaving , ' yes ') ~= 1)title ( ' red T_ {2 r} ; blue T_2 ' ,...' FontName ', police , ' FontSize ', size , ' ← (cid:45) FontWeight ',' bold '); else print ( hFig3 , '- dpdf ', ' ../ GraphicsImages / ← (cid:45) twoLinkArmCtrl2 . pdf ');endhFig4 = figure (4) ; set ( gca , ' FontName ', police , ' FontSize ', size );plot (tr , hx - hxr , 'b ', ' LineWidth ', lineWidth ); grid ;xlabel ( ' time (s) ', ' FontName ', police , ' FontSize ', ← (cid:45) size );ylabel ( ' h_x - h_ { xr } (m) ', ' FontName ', police , ' ← (cid:45) FontSize ', size );if ( strcmp (S. forSaving , ' yes ') ~= 1) title ( ' h_x - h_ { xr } ' ,...' FontName ', police , ' FontSize ', size , ' ← (cid:45) FontWeight ',' bold ');else print ( hFig4 , '- dpdf ', ' ../ GraphicsImages / ← (cid:45) twoLinkArmErrHx . pdf '); end hFig5 = figure (5) ;set ( gca , ' FontName ', police , ' FontSize ', size );plot (tr , hy - hyr , 'b ', ' LineWidth ', lineWidth ); grid ;xlabel ( ' time (s) ', ' FontName ', police , ' FontSize ', ← (cid:45) size ); ylabel ( ' h_y - h_ { yr } (m) ', ' FontName ', police , ' ← (cid:45) FontSize ', size );if ( strcmp (S. forSaving , ' yes ') ~= 1)title ( ' h_y - h_ { yr } ' ,...' FontName ', police , ' FontSize ', size , ' ← (cid:45) FontWeight ',' bold ');else print ( hFig5 , '- dpdf ', ' ../ GraphicsImages / ← (cid:45) twoLinkArmErrHy . pdf ');end% forward reference and actual kinematicshFig6 = figure (6) ; set ( gca , ' FontName ', police , ' FontSize ', size );plot (hx , hy , 'b - ', hxr , hyr , 'r -- ', ' LineWidth ', ← (cid:45) lineWidth ); grid ;xlabel ( ' h_x (m) ', ' FontName ', police , ' FontSize ', ← (cid:45) size );ylabel ( ' h_y (m) ', ' FontName ', police , ' FontSize ', ← (cid:45) size );if ( strcmp (S. forSaving , ' yes ') ~= 1) title ( ' blue h_x h_y red h_ { xr } h_ { yr } ' ,...' FontName ', police , ' FontSize ', size , ' ← (cid:45) FontWeight ',' bold ');else print ( hFig6 , '- dpdf ', ' ../ GraphicsImages / ← (cid:45) twoLinkArmHxHy . pdf ');end hFig7 = figure (7) ;set ( gca , ' FontName ', police , ' FontSize ', size );R. hxi = 0.8*( P. l1 +P. l2 ); R. hxf = 0;R. hyi = P. l1 + 0.1* P. l2 ; R. hyf = - 0.1* P. l1 ; plot ([] , []) ; grid ;hold onfor i = 1:3: length ( ts )line ([0 , l1 * cos ( theta1 (i))], [0 , l1 * sin ( theta1 (i) ← (cid:45) )], ' Color ', [0 0 0.5] , ' LineWidth ', 2) ;line ([ l1 * cos ( theta1 (i)) , l1 * cos ( theta1 (i))+ l2 * cos ← (cid:45) ( theta1 (i)+ theta2 (i)) ] ,... [ l1 * sin ( theta1 (i)) ,... . l1 * sin ( theta1 (i))+ l2 * sin ( theta1 (i)+ theta2 (i ← (cid:45) ))], ' Color ', [0 0 0.9] , ' LineWidth ', ← (cid:45)
2) ;plot ( hx (i) , hy (i) , 'ro ', ' LineWidth ', 2) ; grid ;pause (0.01) ;end hold offxlabel ( ' h_x (m) ', ' FontName ', police , ' FontSize ', ← (cid:45) size );ylabel ( ' h_y (m) ', ' FontName ', police , ' FontSize ', ← (cid:45) size );if ( strcmp (S. forSaving , ' yes ') ~= 1)title ( ' blue h_x h_y red h_ { xr } h_ { yr } ' ,... ' FontName ', police , ' FontSize ', size , ' ← (cid:45) FontWeight ',' bold ');else print ( hFig7 , '- dpdf ', ' ../ GraphicsImages / ← (cid:45) twoLinkArmAnimation . pdf ');end366