Achieving Precise Mechanical Control in Intrinsically Noisy Systems

How can precise control be realised in intrinsically noisy systems? Here, we develop a general theoretical framework that provides a way to achieve precise control in signal-dependent noisy environments. When the control signal has Poisson or supra-Poisson noise, precise control is not possible. If, however, the control signal has sub-Poisson noise, then precise control is possible. For this case, the precise control solution is not a function, but a rapidly varying random process that must be averaged with respect to a governing probability density functional. Our theoretical approach is applied to the control of straight-trajectory arm movement. Sub-Poisson noise in the control signal is shown to be capable of leading to precise control. Intriguingly, the control signal for this system has a natural counterpart, namely the bursting pulses of neurons --trains of Dirac-delta functions-- in biological systems to achieve precise control performance.


Introduction
Many mechanical and biological systems are controlled by signals which contain noise. This poses a problem. The noise apparently corrupts the control signal, thereby preventing precise control. However, precise control can be realised, despite the occurrence of noise, as has been demonstrated experimentally in biological systems. For example, in neural-motor control, as reported in [1], the movement error is believed to be mainly due to inaccuracies of the neural-sensor system, and not associated with the neural-motor system.
The minimum-variance principle proposed in [2,3] has greatly influenced the theoretical study of biological computation. Assuming the magnitude of the noise in a system depends strongly on the magnitude of the signal, the conclusion of [2,3] is that a biological system is controlled by minimising the execution error.
A key feature of the control signal in a biological system is that biological computation often only takes on a finite number of values. For example, 'bursting' neuronal pulses in the neural-motor system control seem very likely to have only three states, namely inactive, excited, and inhibited. This kind of signal (neuronal pulses) can be abstracted as a dynamic trajectory which is zero for most of the time, but intermittently takes a very large value. Generally, this kind of signal looks like a train of irregularly spaced Dirac-delta functions. In this work we shall theoretically investigate the way signals in realistic biological systems are associated with precise control performance. We shall use bursting neuronal pulse trains as a prototypical example of this phenomenon. In a biological system, noise is believed to be inevitable and essential; it is a part of a biological signal and, for example, the magnitude of the noise typically depends strongly on the magnitude of the signal [2,3]. One characteristic of the noise in a system is the dispersion index, α, which describes the statistical regularity of the control signal. When the variance in the control signal is proportional to the 2α-th power of the mean control signal, the dispersion index of the control noise is said to be α. It was shown in [2,3] and elsewhere (e.g., [4,5]) that an optimal solution of analytic form can be found when the stochastic control signal is supra-Poisson, i.e., when α ≥ 0.5. However, the resulting control is not precise and a non-zero execution error arises. In recent papers, a novel approach was proposed to find the optimal solution for control of a neural membrane [6], and a model of saccadic eye movement [7]. It was shown that if the noise of the control signal is more regular than Poisson process (i.e., if it is sub-Poisson, with α < 0.5), then the execution error can be shown to reduce towards zero [6,7]. This work employed the theory of Young measures [13,14], and involved a very specific sort of solution (a 'relaxed optimal parameterized measure solution'). We note that many biological signals are more regular than a Poisson process: e.g., within invivo experiments, it has often been observed that neuronal pulse signals are sub-Poisson in character (α < 0.5) [15,16]. However, in [6,7], only a one-dimensional linear model was studied in detail. Thus the results and methods cannot be applied to the control of general dynamical systems. The work of [6,7] however, leads to a much harder problem: the general mathematical link between the regularity of the signal's noise and the control performance that can be achieved.
In the present work we establish some general mathematical principles linking the regularity of the noise in a control signal with the precision of the resulting control performance, for general nonlinear dynamical systems of high dimension. We establish a general theoretical framework that yields precise control from a noisy controller using modern mathematical tools. The control signal is formulated as a Gaussian (random) process with a signal-dependent variance. Our results show that if the control signal is more regular than a Poisson process (i.e., if α < 0.5), then the control optimisation problem naturally involves solutions with a specific singular character (parameterized measure optimal solutions), which can achieve precise control performance. In other words, we show how to achieve results where the variance in control performance can be made arbitrarily small. This is in clear contrast to the situation where the control signals are Poisson or more random than Poisson (α ≥ 0.5), where the optimal control signal is an ordinary function, not a parameterized measure, and the variance in control performance does not approach zero. The new results can be applied to a large class of control problems in nonlinear dynamical systems of high dimension. We shall illustrate the new sort of solutions with an example of neural-motor control, given by the control of straight-trajectory arm movements, where neural pulses act as the control signals. We show how pulse trains may be realised in nature which lead towards the optimisation of control performance.

Model and Mathematical Formulation
To establish a theoretical approach to the problem of noisy control, we shall consider the where: t is time (t ≥ 0), x(t) = [x 1 (t), . . . , x n (t)] ⊤ is a column vector of 'coordinates' describing the state of the system to be controlled (a ⊤-superscript denotes transpose) and u(t) = [u 1 (t), . . . , u m ] ⊤ , is a column vector of the signals used to control the x system. The dynamical behaviour of the x system, in the absence of a control signal, is determined by a(x, t) and b(x, t), where a(x, t) consists of n functions: a(x, t) = [a 1 (x, t), . . . , a n (x, t)] ⊤ and b(x, t) is an n×m 'gain matrix' with elements b ij (x, t). The system (1) is a generalisation of the dynamical systems studied in the literature [2,3,6,7]. As stated above, the control signal, u(t), contains noise. We follow Harris's work [2, 3] on signal-dependent noise theory by modelling the components of the control signal as where λ i (t) is the mean control signal at time t of the i'th component of u(t) and all noise (randomness) is contained in ζ i (t). In particular, we take the ζ i (t) to be an independent Gaussian white noises obeying is Dirac-delta function, and δ ij is Kronecker delta. The quantities σ i (t), which play the role of standard deviations of the ζ i (t), are taken to explicitly depend on the mean magnitudes of the control signals: where κ i is a positive constant and α is the dispersion index of the control process (described above). Thus, we can formulate the dynamical system, Eq. (1), as a system of Itô diffusion equations: , · · · , n and j = 1, . . . , m. We make the assumption that the range of each ξ i is bounded: be the region where the control signal takes values, with m denoting the m-order Cartesian product. Let Ξ be state space of x. In this paper we assume it be bounded.
Let us now introduce the function φ(x, t) = [φ 1 (x, t), · · · , φ k (x, t)] ⊤ , which represents the objective that is to be controlled and optimised. For example, for a linear output we can take φ(x, t) = Cx for some k ×n matrix C; in the case that we control the magnitude of x, we can take φ(x, t) = x 2 ; we may even allow dependence on time, for example if the output decays exponentially with time, we can take exp(−γt)x(t) for some constant γ > 0.
The aim of the control problem we consider here is: (i) to ensure the expected trajectory of the objective φ(x(t), t) reaches a specified target at a given time, T , and (ii) to minimise the execution error accumulated by the system, during the time, R, that the system is required to spend at the 'target' [2,3,6,7,8,9,10,11]. In the present context, we take the motion to start at time t = 0 and subject to the initial condition x(0). The target has coordinates z(t) and we need to choose the controller, u(t), so that for the time interval T ≤ t ≤ T + R the expected state of the objective φ(x(t), t) of the system satisfies E[φ(x(t), t)] = z(t). The accumulated execution error is T +R T i var (φ i (x(t), t)) dt and we require this to be minimised. Statistical properties of x(t) can be written in terms of p(x, t), the probability density of the state of the system (4) at time t, which satisfies the Fokker-Plank equation with Three important quantities are the following: (A) The accumulated execution error: (C) The dynamical equation of p(x, t) described as (5).

The Young Measure Optimal Solution
To illustrate the idea of the solutions we introduce here, namely Young measure optimal solutions, we provide a simple example. Consider the situation where x and u are onedimensional functions, while a(x, t) = px, b(x, t) = q, κ = 1, z(t) = z 0 and φ(x, t) = x.
This has the solution )qλ(s)ds and its variance is var(x(t)) = t 0 exp(2p(t − s))q 2 |λ(s)| 2α ds. The solution of the optimisation problem is the minimum of the following functional: with for some integrable function γ(t), which serves as a Lagrange auxiliary multiplier.
In the general case, we minimise (A) using (B) and (C) as constraints via the introduction of appropriate x and t dependent Lagrange multipliers. This leads to a functional of the mean control signal, H [λ], with the form H[λ] = T +R T h(t, λ(t))dt (see below and Appendix A).
Let us use ξ = [ξ 1 , · · · , ξ m ] ⊤ to denote the value of λ(t) at a given time of t, i.e., ξ = λ(t); ξ will serve as a variable of the Young measure (see below). We find where g i (t), f i (t) and z(t) are functions with respect to t but are independent of the variable ξ.
The abstract Hamiltonian minimum (maximum) principle (AHMP) [12] provides a necessary condition for the optimal solution of minimising (A) with (B) and (C), which is composed of the points in the domain of definition of λ, namely, Ω, that minimize the function h(t, ξ) in (8), at each time, t, which is named Hamiltonian integrand. This principle tells us that the optimal solution should pick values of the minimum of h(t, ξ) with respect to ξ, for each t.
If the control signal is supra-Poisson or Poisson, namely the dispersion index α ≥ 0.5, for each t ∈ [0, T + R], the Hamiltonian integrand h(t, ξ) is convex (or semi-convex) with respect to ξ and so has a unique minimum point with respect to each ξ i . So, the optimal solution is a deterministic function of time: for each t 0 , λ i (t 0 ) can be regarded as picking value at the minimum point of h(t, ξ) for t = t 0 .
When α < 1/2, namely when the control signal is sub-Poisson, it follows that h(t, ξ) is no longer a convex function. Figs. 1 show the possible minimum points of the term g i (t)|ξ i | 2α − f i (t)ξ i with g i (t) > 0 and f i (t) > 0. From the assumption that the range of each ξ i is bounded, namely −M Y ≤ ξ i ≤ M Y , it then directly follows, from the form of h(t, ξ), that the value of ξ i which optimises h(t, ξ) is not unique; there are three possible minimum values: −M Y , 0, and M Y , as shown in Table 1. So, no explicit function λ(t) exists which is the optimal solution of the optimisation problem (A)-(C). However, an infinimum of (8) does exist.
Proceeding intuitively, we first make an arbitrary choice of one of the three optimal values for ξ i (namely one of −M Y , 0, and M Y ) and then average over all possible choices at each time. With η t,i (ξ) the probability density of ξ i at time t, the average is carried out using the distribution (probability density functional) η[λ] ∝ t,i η t,i (ξ) which represents independent choices of the control signal at each time. Thus, for example, the functional H[λ] becomes functionally averaged over λ(·) according to The optimisation problem has thus shifted from determining a function (as required when α ≥ 1/2) to determining a probability density functional, η[λ]. This intuitively motivated procedure is confirmed by optimisation theory-and this leads us to Young measure theory.
Let us spell it out in a mathematical way. Young measure theory [13,14] provides a solution to an optimization problem where a solution, which was a function, becomes a linear functional of a parameterized measure. By way of explanation, a function, λ(t), yields a single value for each t, but a parameterized measure {η t (·)} yields a set of values on which a measure (i.e., a weighting) η t (·) is defined for each t. A functional with respect to a parameterized measure can be treated in a similar way to a solution that is an explicit function, by averaging over the set of values of the parameterized measure at each t. In detail, a functional of the form H[λ] = T 0 h(t, λ(t))dt, of an explicit function, λ(t), can have its definition extended to a parameterized measure η t (·), namely H In this sense, an explicit function can be regarded as a special solution that is a 'parameterized concentrated measure' (i.e., involving a Dirac-delta function) in that we can write Thus, we can make the equivalence between the explicit function λ(t) and a parameterized concentrated measure {δ(ξ − λ(t))} t and then replace this concentrated measure, when appropriate, by a Young measure.
Technically, a Young measure is a class of parameterized measures that are relatively weak*-compact such that the Lebesgue function space can be regarded as its dense subset in the way mentioned above. Thus, by enlarging the solution space from the function space to the (larger) Young measure space, we can find a solution in the larger space and the minimum value of the optimisation problem, in the Young measure space, coincides with the infinimum in the Lebesgue function space.
For any function r(x, t, ξ), we denote a symbol · as the inner product of r(x, t, ξ) over the parameterized measure η t (dξ), by averaging r(x, t, ξ) with respect to ξ via η t (·). That is we define r(x, t, ξ)·η t to represent Ω r(x, t, ξ)η t (dξ). In this way we can rewrite the optimisation problem (A)-(C) as: Here, Y denotes the Young measure space, which is defined on the state space Ξ with t ∈ [0, T + R], while η = {η t (·)} denotes a shorthand for the Young measure associated with control; (L · η) • p is defined as So, we can study the relaxation problem (9) instead of the original one, (A)-(C). We assume that the constraints in (9) admit a nonempty set of λ(t), which guarantees that the problem (9) has a solution. We also assume the existence and uniqueness of the Cauchy problem of the Fokker-Plank equation (5). The abstract Hamiltonian minimum (maximum) principle (Theorem 4.1.17 [12]) also provides a similar necessary condition for the Young measure solution of (9), if it admits a solution, that is composed of the points in Ω which minimise the integrand of the underlying 'abstract Hamiltonian'. By employing variational calculus with respect to the Young measure, we can derive the form (8), for the Hamiltonian integrand. See Appendix A for details.
Via this principle, the problem conceptively reduces to finding the minimum points of h(t, ξ). From Table 1, for a sufficiently large M Y , it can be seen that, if α < 0.5, then the minimum points for each t with Hence, in the case of α < 0.5, the optimal solution of (9) is a measure on {M Y , 0} or {−M Y , 0}. This implies that the optimal solution of (9) should have the following form η t (·) = η 1,t (·)×, · · · , η m,t (·), where × stands for the Cartesian product, and each η i,t we adopt where µ i (t) and ν i (t) are non-negative weight functions. The optimisation problem corresponds to the determination of the µ i (t) and ν i (t). Averaging with respect to η corresponds to the optimal control signal when the noise is sub-Poisson (α < 0.5). This assignment of a probability density for the solution at each time is known in the mathematical literature as a Young Measure [12,13,14]. For all i and t, the weight functions satisfy: . Consider the simple one-dimensional system (6). We shall provide the explicit form of the optimal control signal u(t) as a Young measure. Taking expectation for both sides in (6), Since we only minimise the variance in [T, T + R] for some T > 0 and R > 0, the control signal u(t) for t ∈ [0, T ) is picked so that the expectation of x(t) can reach z 0 at the time t = T . After some simple calculations, we find a deterministic λ(t) as follows: In the interval [T, T + R], as discussed above, for a sufficiently large M Y , the optimal solution of λ(t) should be a Young measure that picks values in {0, M Y , −M Y }. To sum up, we can construct the optimal λ(t) as follows: It can be seen that in [0, T ), η t (·) is in fact a deterministic function as the same as λ(t).

Precise Control Performance
We now illustrate the control performance when the noise is sub-Poisson. For the general nonlinear system (1), we cannot obtain an explicit expression for the probability density functional η[λ], Eq. (10), or the value of the variance (execution error). However, we can adopt a non optimal probability density functional which illustrates the property of the exact system, that the execution error becomes arbitrarily small when the bound of the control signal, M Y , becomes arbitrarily large. In the simple case (6), we note that if there is aû(t), such that E(x(t)) = z(t), then the variance becomes, expressed by Young measureη(·), which converges to zero as M Y → ∞, due to α < 0.5. That is, the minimised execution error can be arbitrarily small if the bound of the control signal, M Y , goes sufficiently large.
In fact, this phenomenon holds for general cases. The non optimal probability density functional is motivated by assuming that there is a deterministic control signalû(t) which controls the dynamical system which is the original system (1), with the noise removed. The deterministic control signalû(t) causesx(t) to precisely achieve the target trajectoryx(t) = z(t) for T ≤ t ≤ T + R. Then, we add the noise with the signal-dependent variance: σ i = κ i |λ i | α with some α < 0.5, which leads a stochastic differential equation, dx = A(x, t, λ(t))dt + B(x, t, λ(t))dW t . The non optimal probability density that is appropriate for time t, namelyη t,i (ξ), is constructed to have a mean over the control values {−M Y , 0, M Y }, which equalsû(t). This probability density iŝ where σ(t) = sign(û i (t)) and, by definition, We establish in Appendix B that the expectation condition ((B) above) holds asymptotically when M Y → ∞, which shows that the non optimal probability density functional is appropriately 'close' to the optimal functional. The accumulated execution error associated with the non optimal functional is estimated as and, in this way, optimal performance of control, with sub-Poisson noise, can be seen to become precise as M Y is made large. By contrast, if α ≥ 0.5, the accumulated execution error is always greater than some positive constant.
To gain an intuitive understanding of why the effects of noise are eliminated for α < 0.5 we discretise the time t into small bins of identical size ∆t. Using the 'noiseless control'û i (t), we divide the time bin [t, t + ∆t] into two complementary intervals: [t, t + |û(t)|∆t/M Y ] and [t + |û(t)|∆t/M Y , t + ∆t], and assign λ i = σ(t)M Y for the first interval and λ i = 0 for the second. When ∆t → 0 the effect of the control signal λ i (t) on the system approaches that ofû i (t), although λ i (t) andû i (t) are quite different. The variance of the noise in the first interval is κ i M 2α Y and is 0 in the second. Hence, the overall noise effect of the bin is Remarkably, this tends to zero as M Y → ∞ if α < 1/2 (i.e., for sub-Poisson noise). The discretisation presented may be regarded as a formal stochastic realisation of the probability density functional (Young measure) adopted. The interpretation above can be verified in a rigorous mathematical way. See Appendix B for details.

Application and Example
Let us now consider an application of this work: the control of straight-trajectory arm movement, which has been widely studied [8,9,10,11] and applied to robotic control. The dynamics of such structures are often formalised in terms of coordinate transformations. Nonlinearity arises from the geometry of the joints. The change in spatial location of the hand that results from bending the elbow depends not only on the amplitude of the elbow movement, but also on the state of the shoulder joint.
For simplicity, we ignore gravity and viscous forces, and only consider the movement of a hand on a horizontal plane in the absence of friction. Let θ 1 denote the angle between the upper arm and horizontal direction, and θ 2 be the angle between the forearm and upper arm (Fig. 2). The relation between the position of hand [x 1 , x 2 ] and the angles [θ 1 , θ 2 ] is θ 1 = arctan(x 2 /x 1 ) − arctan(l 2 sin θ 2 /(l 1 + l 2 cos θ 2 )) , where l 1,2 are moments of inertia with respect to the center of mass, for the upper arm and forearm. When moving a hand between two points, a human maneuvers their arm so as to make the hand move in roughly a straight line between the end points. We use this to motivate the model by applying geostatics theory [8]. This implies that the arm satisfies an Euler-Lagrange equation, which can be described as the following nonlinear two-dimensional system of differential equations: In these equations N =    I 1 + m 1 r 2 1 + m 2 l 2 1 +I 2 + m 2 r 2 2 + 2k cos θ 2 I 2 + m 2 r 2 2 + k cos θ 2 where m i , l i , and I i are, respectively the mass, length, and moment of inertia with respect to the center of mass for the i'th part of the system and i = 1 (i = 2) denotes the upper arm (forearm), r 1,2 are the lengths of the upper-and fore-arms, and γ 0 is the scale parameter of the force. Additionally, k = m 2 l 1 r 2 , while λ 1,2 (t) are the means of two torques Q 1,2 (t), which are motor commands to the joints. The torques are accompanied by signal-dependent noises. All other quantities are fixed parameters. See [8] for the full details of the model. The values of the parameters we pick here are listed in Table 2.
For this example, we shall aim to control the hand such that it starts at t = 0, with the initial condition of (14), reaches the target at coordinates H = [H 1 , H 2 ] at time t = T , and then stays at this target for a time interval of R. We use the minimum variance principle to determine the optimal task, which is more advantageous than other optimisation criteria to control a robot arm [8,11]. Let [x 1 (t), x 2 (t)] be the Cartesian coordinates of the hand that follow from the angles [θ 1 (t), θ 2 (t)]. The minimum variance principle determines min λ 1 ,λ 2 Despite not being in possession of an explicit analytic solution, we can conclude that if α ≥ 0.5, the optimisation problem results from the unique minimum to the Hamiltonian integrand and hence yields λ 1 (t) and λ 2 (t) which are ordinary functions. However, if α < 0.5, the optimal solution of the optimisation problem follows from a probability density functional analogous to Eq. (10) (i.e., a Young measure over λ i ∈ {−M Y , 0, M Y }). Thus, we can relax the optimisation problem via Young measure as follows: We used Euler's method to conduct numerical computations, with a time step of 0.01 msec in (14). This yields a dynamic programming problem (see Methods). Fig. 3 shows the means of the optimal control signalsλ 1,2 (t) with α = 0.25 and M Y = 20000: According to the form of the optimal Young measure, the optimal solution should be It can be shown (derivation not given in this work) that in the absence of the noise term, the arm can be accurately controlled to reach a given target for any T > 0. In this case, Fig. 4 shows the dynamics of the angles, their velocities, and accelerators, in the controlled system, removed noise. See, in comparison, the dynamical system with noise, whose dynamics of the angles, velocities, and accelerations are illustrated in Fig. 5, and its dynamics are exactly the same as those in the case with noise removed. However, the acceleration dynamics of a noisy dynamic system appear discontinuous since the control signals, that have noises and are added to the right-hand sides of the mechanical equations, are discontinuous (noisy) in a numerical realisation. However, according to the theory of stochastic differential equations [17], (14) has continuous solution. Hence, these discontinuous acceleration dynamics lead very smooth dynamics of velocities and angles, as shown in Fig. 5.
Figs. 6 (a) and (b) illustrate that the probability density functional, for this problem, contains optimal control signals that are similar to neural pulses. Despite the optimal solution not being an ordinary function when α < 0.5, the trajectories of the angles θ 1 and θ 2 of the arm appear quite smooth, as shown in Fig. 5 (a), and the target is reached very precisely if the value of M Y is large. By comparison, when α > 0.5 the outcome has a standard deviation between 4 to 6 cm, which may lead to a failure to reach the target. A direct comparison between the execution error of the cases α = 0.8(> 0.5) and α = 0.25(< 0.5) is shown in the supplementary movies (supplementary videos 'Video S1' and 'Video S2') of arm movements of both cases. Our conclusion is that a Young measure optimal solution, in the case of sub-Poisson control signals, can realize a precise control performance even in the presence of noise. However, Poisson or Supra-Poisson control signals cannot realise a precise control performance, despite the existence of an explicit optimal solution in this case. Thus α < 0.5 significantly reduces execution error compared with α ≥ 0.5.
With different T (the starting time of reaching the target) and R (the duration of reaching the target), under sub-Poisson noise, i.e., α < 0.5, the system can be precisely controlled by optimal Young measure signals with a sufficiently large M Y . Since the target in the reachable region of the arm, it implies that the original differential system of (14) with the noise removed can be controlled for any T > 0 and R > 0 [8,9]. According to the discussion in Appendix B (Theorem 2), the execution error can be arbitrarily small when M Y is sufficiently large. However, for a smaller T , i.e., the more rapid the control is, the larger means of the control signals will be. As for the duration R, by picking the control signals as fixed values (zeros in this example) such that the velocities keep zeros, the arm will stay at the target for arbitrarily long or short. Similarly, with a large M Y , the error (variance) of staying at the target can be very small. To illustrate these arguments, we take T = 100 (msec) and R = 100 (msec) for example (all other parameters are the same as above). Fig. 7 shows that the means of the optimal Young measure control signals before reaching the target have larger amplitudes than those when T = 650 (msec) and Fig. 8 shows that the arm can be precisely controlled to reach and stay at the target.
The movement error depends strongly on the value of the dispersion index, α, and the bound of the control signal, M Y . Fig. 9 indicates a quantitative difference in the execution error between the two cases α < 0.5 and α ≥ 0.5, if α is close to (but less than) 0.5. The execution error can be appreciable unless a large M Y is used. For example if α = 0.45, as in Fig. 9, the square root of the execution error is approximately 0.6 cm when M Y = 20000. From (13), the error decreases as M Y increases, behaving approximately as a power-law, as illustrated in the inner plot of Fig. 9. The logarithm of the square root of the execution error is found to depend approximately linearly on the logarithm of M Y when α = 0.25, with a slope close to −0.25, in good agreement with the theoretical estimate (13).
We note that in a biological context, a set of neuronal pulse trains can achieve precise control in the presence of noise. This could be a natural way to approximately implement the probability density functional when α < 0.5. All other parameters are the same as above (α = 0.25). The firing rates are illustrated in Fig. 6 (a) and (b) and broadly coincide with the probability density functional we have discussed. In particular, at each time t, the probability η i,t can be approximated by the fraction of the neurons that are firing, with the mean firing rates equal the means of the control signals (see Methods). The approximations of the components of the noisy control signals are shown in Figs. 10 (a) and (b) respectively. Fig.10 (c) and (d) illustrate such an implementation of the optimal solution by neuronal pulse trains. Using the pulse trains as control signals, we can realise precise movement control. We enclose two videos 'movieUP.avi' and 'movieDOWN.avi' to demonstrate the efficiency of the control by pulse trains with two different targets. As they show, the targets are precisely accessed by the arm. We point out that the larger the ensemble is, the more precise the control performance will be, because a large number of the neurons in an ensemble can theoretically lead to a large M Y as we mentioned above, which results in an improvement of the approximation of a Young measure and decreases the execution error as stated in (13).
We note that these kinds of patterns of pulse trains have been widely reported in experiments, for example, the synchronous neural bursting reported in [18]. This may provide a mathematical rationale for the nervous system to adopt pulse-like signals to realise motor control.

Conclusions
In this paper, we have provided a general mathematical framework for controlling a class of stochastic dynamical systems with random control signals whose noisy variance can be regarded as a function of the signal magnitude. If the dispersion index, α, is < 0.5, which is the case when the control signal is sub-Poisson, an optimal solution of explicit function does not exist but has to be replaced by a Young measure solution. This parameterized measure can lead a precise control performance, where the controlling error can become arbitrarily small. We have illustrated this theoretical result via a widely-studied problem of arm movement control.
In the control problem of biological and robotic systems, large control signals are needed for rapid movement control [21]. When noise occurs, this will cause imprecision in the control performance. As pointed out in [9,22,23,24,25], a trade-off should be considered when conducting rapid control with noises. In this paper, we still use a "large" control signal but with different contexts. With sub-Poisson noises, we proved that a sufficiently large M Y , i.e., a sufficiently large region of the control signal values, can lead precise control performance. Hence, a large region of control signal values plays a crucial role in realising precise control in noisy environments, for both "slow" and "rapid" movement control. In numerical examples, the larger M Y we pick, the smaller control error will be, as shown in the inset plot of Fig. 9 as well as (13) (Theorem 2 in Appendix B).
Implementation of the Young measure approach in biological control appears to be a natural way to achieve precise execution error in the presence of sub-Poisson noise. In particular, in the neural-motor control example illustrated above, the optimal solution in the case of α < 0.5 is quite interesting. Assume we have an ensemble of neurons which fire pulses synchronously within a sequence of non overlapping time windows, as depicted in Figs. 6 C and D. We see that the firing neurons yield control signals which are very close, in form, to the type of Young measure solution. This conclusion may provide a mathematical rationale for the nervous system why to adopt pulse-like trains to realise motor control. Additionally, we point out that, our approach may have significant ramifications in other fields, including robot motor control and sparse functional estimation, which are issues of our future research.

Numerical methods for the optimisation solution.
We used Euler's method to conduct numerical computations, with a time step of ∆t = 0.01 msec in (14) with α < 0.5. This yields a dynamic programming problem. First, we divide the time domain [0, T ] into small time bins with a small size ∆t. Then, we regard the process η i,t in each time bin [n∆t, (n + 1)∆t] as a static measure variable. Thus, the solution reduces to finding two series of nonnegative parameters µ i,n and ν i,n with µ i,n + ν i,n ≤ 1 and µ i,n ν i,n = 0 such that n∆t ≤ t < (n + µ i,n )∆t −M Y (n + µ i,n )∆t ≤ t < (n + µ i,n + ν i,n )∆t 0 (n + µ i,n + ν i,n )∆t ≤ t < (n + 1)∆t. .
The approximate solution of the optimisation problem requires nonnegative µ i,n and ν i,n that minimise the final movement errors. We thus have a dynamic programming problem. We should point out that in the literature, a similar method was proposed to solve the optimisation problem in a discrete system with control signals taking only two values [11,19,20]. Thus, the dynamical system (1) becomes the following difference equations via the Euler method: where ν j , j = 1, . . . , m, are independent standard Gaussian random variables. We can derive difference equations for the expectations and variances of x(k), by ignoring the higher order terms with respect to ∆t: Thus, Eq. (9) becomes the following discrete optimization problem: with x(k) a Gaussian random vector with expectation E(x(k)) and covariance matrix cov(x(k), x(k)). Neuronal pulse trains approximating Young measure solution. At each time t, the measure η * t can be approximated by the fraction of the neuron ensemble that are firing. In detail, assuming that the means of the optimal control signals are R t,i (ξ), i = 1, 2, and there are one ensemble of excitatory neurons and another ensemble of inhibitory neurons. A fraction of the neurons fire so that the mean firing rates satisfy: where R ext i (t) and R inh i (t) are the firing rates of the excitatory and inhibitory neurons respectively, and γ is a scalar factor. In occurrence of sub-Poisson noise, the noisy control signals u * . Both ensembles of neurons are imposed with baseline activities, which bound the minimum firing rates away from zeros, given the spontaneous activities of neurons when no explicit signal is transferred. A numerical approach involves discretise time t into small bins of identical size ∆t. The firing rates can be easily estimated by averaging the population activities in a time bin. We have used 400 neurons to control the system, with two ensembles of neurons with equivalent numbers that approximate the first and second components of control signal, respectively. Each neuron ensemble have 200 neurons with 100 excitatory and 100 inhibitory neurons.

Appendices
Appendix A: Derivation of formula (8) Let: W be the time-varying p.d.f. p(x, t) that is second-order continuous-differentiable with respect to x and t that is embedded in the Sobolev function space W 2,2 ; W 1 be the function space of p(x, t 0 ), regarding as a function with respect to x with a fixed t 0 ; W 2 be the function space of p(x 0 , t), regarding as a function with respect to t with a given x 0 . The spaces W 1,2 can be regarded being embedded in W . In addition, let:Ŵ be the function space where the image L[W ] is embedded; L be the space of linear operator L, denoted above; L (Z, E) be the space composed of bounded linear operator from linear space Z to E; and Z * be the dual space of the linear space Z: Z * = L (Z, R). Furthermore, letỸ be the tangent space of Young measure space Y:Ỹ = {η − η ′ : η, η ′ ∈ Y}. For simplicity, we do not specify the spaces and just provide the formalistic algebras, and then the following is similar to Chapter 4.3 in [12] with appropriate modifications. Define Thus, (9) can be rewritten as: The Gâteaux differentials of these maps with respect to p(x, t), denoted by ∇ p ·, are: . And, the differentials of these maps with respect to the Young measure η are: for two Young measuresη, η ∈ Y. Here, ∇ η Φ ∈Ỹ * , ∇ η Π ∈ L (Ỹ,Ŵ × W 1 ), and ∇ η J ∈ L (Ỹ, W 2 ). Then, we are in the position to derive the result of (8) by the following theorem, as a consequence from Theorem 4.1.17 in [12]. (17). Assume that: (1). the trajectory of x(t) in (4) is bounded almost surely; (2). a(x, t), b(x, t) and φ(x, t) are C 2 with respect to (x, t). Let (η * , p * ) be the optimal solution of (9). Then, there are some λ 1 ∈ L (Ŵ × W 1 , W * ), λ 2 = [λ 21 , λ 22 ] ⊤ with λ 21 ∈ L (Ŵ , W 2 ), and λ 22 ∈ L (W 1 , , W 2 ), such that and the abstract maximum principle holds with "abstract Hamiltonian": Proof. Under the conditions in this theorem, we can conclude that the Fokker-Planck equation has a unique solution p(η) that is continuously dependent of η from theory of stochastic differential equation [17]; Π(·, η) : W → W * is Fréchet differentiable at p = π(η) because x(t) is assumed to be almost surely bounded; Π(p, ·) : Y → W * and J (in fact [26]; the partial differential ∇ η Π is weak-continuous with respect to η because it is linearly dependent of η. In addition, from the existence and uniqueness of the Fokker-Planck equation, ∇ p Π(p, η) : W → Im(∇ p Π(p, η)) ⊂W × W 1 has a bounded inverse. This implies that the followingadjoint equation has a solution for p = p(η), denoted by µ. Let µ = [µ 1 , µ 2 ], which should be solution of the following equation with the dual operator L * of L (the operator in the back-forward Kolmogorov equation), still dependent of (x, t) and the value of λ(t) (namely ξ in Young measure): We pick λ 1,2 (t) with µ 1 = λ 1 +c(t)λ 21 and µ 2 = λ 1 +c(t)λ 22 . So, λ 1,2 should satisfy equation (19). In fact, λ 1,2 can be regarded as functions ( or generalized functions) with respect to (x, t).
Thus, the conditions of Lemma 1.3.16 in [12] can be verified, which implies that the gradients of the maps Φ and J with respect to η, by regarding p = p(η) from Π(p, η) = 0, as follows: From the abstract Hamilton minimum principle (Theorem 4.1.17 in [12]), applied to each solution of (9), denoted by η * , there exists a nonzero function c(t) such that is an 'abstract Hamiltonian', with respect toη. With the definitions of λ 1,2 , (23) becomes By specifying µ with λ 1,2 , we have where p * stands for the time-varying density corresponding to the optimal Young measure solution η * . From this, letting ξ = λ, we have the "abstract Hamiltonian" in the form of (21) as the Hamiltonian integrand of H(·). The Hamiltonian abstract minimum (maximum) principle indicates the optimal Young measure η * t is only concentrated at the minimum points of h(t, ξ) with respect to ξ for each t, namely. That is, (20) holds. This completes the proof.
From this theorem, since the variances depend on the magnitude of the signal as described in (3), removing the terms without ξ, it is equivalent to look at the minima of h(t, ξ) in the form of instead of h(t, ξ), where This gives formula (8).
Appendix B: Derivation of precise control performance (13) The control performance inequality (13) can be derived from the following theorem.
In addition, also approaches zero as M Y goes to infinity. This proves the first item of the theorem. This completes the proof.
Hence, as M Y goes to infinity, the non optimal solution (12) can asymptotically satisfy the constraint and the error variance goes to zero as M Y goes to infinity. Therefore, the performance error of the REAL optimal solution of the optimisation problem (9) approaches zero as M Y → ∞ in the case of α < 0.5. Furthermore, we can conclude from (25) that the execution error, measured by the standard deviation, can be approximated as (13). where P is fixed and others not, and two arms (upper arm P Q and the forearm QH). Button H is to reach some given target (red cross) by moving front-and back-arms. Table 1. Summary of the possible minimum points of (8).      Table 2. Target is set by θ 1 (T ) = −1 and θ 2 (T ) = π 2 but without noise: The dynamics of the angles (a), the angle velocities (b) and accelerations (c) (the blue solid curves for those of θ 1 and the green solid curves for θ 2 ). The blue and red dash vertical lines stand for the start and end time points of the duration of reaching the target.

Parameters
Values masses (of the inertia w.r.t the mass center) m 1 = 2.28 kg, m 2 = 1.31 kg lengths (of the inertia w.r.t the mass center) l 1 = 0.305 m, l 2 = 0.254 m moments (of the inertia w.r.t the mass center) I 1 = 0.022 kg · m 2 ; I 2 = 0.0077 kg · m 2 lengths of arms    7. Means of the optimal control signals λ 1 (t) (blue solid) and λ 2 (t) (green solid) in the straight-trajectory arm movement example with T = 100 (sec) and R = 100 (msec). The blue and red dash vertical lines stand for the start and end time points of the duration of reaching the target.