Binary Control and Digital-to-Analog Conversion Using Composite NUV Priors and Iterative Gaussian Message Passing
BBINARY CONTROL AND DIGITAL-TO-ANALOG CONVERSION USING COMPOSITE NUVPRIORS AND ITERATIVE GAUSSIAN MESSAGE PASSING
Raphael Keusch, Hampus Malmberg, and Hans-Andrea Loeliger
ETH Zurich, Dept. of Information Technology & Electrical Engineering { keusch, malmberg, loeliger } @isi.ee.ethz.ch ABSTRACT
The paper proposes a new method to determine a binary con-trol signal for an analog linear system such that the state, orsome output, of the system follows a given target trajectory.The method can also be used for digital-to-analog conversion.The heart of the proposed method is a new binary-enforcing NUV prior (normal with unknown variance). Theresulting computations, for each planning period, amount toiterating forward-backward Gaussian message passing re-cursions (similar to Kalman smoothing), with a complexity(per iteration) that is linear in the planning horizon. In conse-quence, the proposed method is not limited to a short planninghorizon.
Index Terms — Discrete-level priors, normals with un-known variance (NUV), finite-control-set model predictivecontrol (MPC), digital-to-analog conversion (DAC).
1. INTRODUCTION
Consider the classical control problem of steering an analogphysical linear system along some desired trajectory, or tomake the system produce some desired analog output. Inthis paper, we are interested in the special case where thecontrol input is binary (i.e., restricted to two levels), whichmakes the problem much harder. This binary-input controlproblem includes, in particular, a certain type of analog-to-digital converter where the binary output of some digital pro-cessor directly drives a continuous-time analog linear filter—preferably an inexpensive one—which produces the desiredanalog waveform.It is tempting to ask for an optimal binary control signal,i.e., a control signal that produces the best approximation ofthe desired analog trajectory (e.g., for a quadratic cost func-tion). However, determining such an optimal control signalis a hard combinatorial optimization problem, with a compu-tational complexity growing exponentially with the planninghorizon [1, 2]. In consequence, insisting on an optimal con-trol signal limits us to a short planning horizon, which is avery severe restriction. This problem is well known in model The method of this paper can be extended to control signals with
M > levels, as will be described elsewhere. predictive control (MPC) [3, 4]. Techniques such as spheredecoding do help [5], but the fundamental problem remains.Clearly, the binary-input control problem is a nonconvexoptimization problem. A general approach to nonconvex opti-mization is to resort to some convex relaxation, and to projectthe solution back to the permissible set [6]. Other generalapproaches include heuristic methods such as random-restarthill-climbing [7] and simulated annealing [8].The heart of the method proposed in this paper is a newbinary-enforcing NUV prior (normal with unknown vari-ance). NUV priors are a central idea of sparse Bayesianlearning [9, 10, 11, 12], and closely related to variationalrepresentations of Lp norms [13, 14]. Such priors have beenused mainly for sparsity; in particular, no binary-enforcingNUV prior seems to have been proposed in the prior liter-ature. (An interesting non-NUV binary-enforcing prior hasbeen proposed in [15].)A main advantage of NUV priors in general is theircomputational compatibility with linear Gaussian models,cf. [16]. In this paper, the computations (for each planningperiod) amount to iterating forward-backward Gaussian mes-sage passing recursions similar to Kalman smoothing, with acomplexity (per iteration) that is linear in the planning hori-zon. In consequence, the proposed method can effectivelyhandle long planning horizons, which can far outweigh itssuboptimality.
2. THE BINARY-ENFORCING NUV PRIOR
Let N (cid:0) x ; µ, σ (cid:1) denote the normal probability density func-tion in x with mean µ ∈ R and variance σ . Let ρ ( x, θ ) (cid:44) N (cid:0) x ; a, σ (cid:1) N (cid:0) x ; b, σ (cid:1) , (1)where θ (cid:44) ( σ , σ ) is a shorthand for the two unknown vari-ances in (1). For fixed θ , ρ ( x, θ ) is a normal probability den-sity in x (up to a scale factor), which is essential for the algo-rithms in Section 3.3.The heart of the proposed method is the observationthat (1) can be used as a (improper) joint prior for x and θ that strongly encourages x to lie in { a, b } . To see this, assume ρ ( x, θ ) is used in some model with (fixed) observation(s) ˘ y © 2021 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating newcollective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. a r X i v : . [ ee ss . SP ] F e b σ × (cid:63) N (cid:45) + (cid:63) a (cid:45)(cid:45) σ × (cid:54) N (cid:45) + (cid:54) b (cid:45) = ρ ( x, θ ) (cid:45) X N (cid:0) , s (cid:1) (cid:63) + (cid:45) µp (˘ y | x ) Fig. 1 . Factor graph of (3) for fixed ˘ y . The boxes labeled “ N ”represent normal probability density functions N (0 , . − . . . µ ˆ x AM s = 0 s = 0 . s = 5 Fig. 2 . The estimate of Section 2.1, for a = 0 and b = 1 , as afunction of µ . − . . . µ ˆ x EM s = 0 s = 0 . s = 5 Fig. 3 . The estimate of Section 2.2, for a = 0 and b = 1 , as afunction of µ .and likelihood function p (˘ y | x ) . Two different ways to esti-mate the variances θ are considered in Sections 2.1 and 2.2.For the numerical examples in Figs. 2 and 3, we will assumethat p (˘ y | x ) is Gaussian (in x ) with mean µ and variance s depending on ˘ y , i.e., p (˘ y | x ) = N (cid:0) x ; µ, s (cid:1) . (2)A factor graph [17] of the resulting statistical system model p (˘ y | x ) ρ ( x, θ ) = N (cid:0) x ; µ, s (cid:1) N (cid:0) x ; a, σ (cid:1) N (cid:0) x ; b, σ (cid:1) (3)is shown in Fig. 1. Assume that x and θ are determined by joint MAP estimation.The resulting estimate ˆ x of x is ˆ x = argmax x max θ p (˘ y | x ) ρ ( x, θ ) (4) = argmax x p (˘ y | x ) max θ ρ ( x, θ ) , (5)with an effective (improper) prior max θ ρ ( x, θ ) = max σ N (cid:0) x ; a, σ (cid:1) max σ N (cid:0) x ; b, σ (cid:1) (6) ∝ | x − a | · | x − b | (7)It is obvious that this effective prior has a strong preferencefor x to lie in { a, b } .Such estimates may conveniently be computed by alter-nating maximization (AM) over x and θ , cf. Section 3.3.2.However, AM may converge to a local maximum rather thanthe global maximum (5), in which case the estimate ˆ x AM re-turned by AM may lie outside { a, b } .For p (˘ y | x ) as in (2), the estimate ˆ x AM is plotted in Fig. 2.In this setting, we observe (and it can be proved) that, forevery fixed µ and sufficiently large s , there is no local maxi-mum and ˆ x AM lies in { a, b } . In this case, we first determine the MAP estimate ˆ θ of θ , i.e., ˆ θ = argmax θ (cid:90) ∞−∞ p (˘ y | x ) ρ ( x, θ ) dx, (8)and then we estimate x as ˆ x = argmax x p (˘ y | x ) ρ ( x, ˆ θ ) . (9)Such estimates may conveniently be computed by expec-tation maximization (EM) (cf. Section 3.3.3), which, how-ever, may converge to a local maximum. Moreover, even theglobal maximum (9) may lie outside { a, b } . Nonetheless, for p (˘ y | x ) as in (2), we observe (and it can be proved) that, forevery fixed µ (cid:54) = ( a + b ) / and sufficiently large s , ˆ x EM liesin { a, b } , as illustrated in Fig. 3.
3. SYSTEM MODEL AND ALGORITHMS3.1. Problem Statement
Consider a linear system with scalar input u k ∈ R and state x k ∈ R N , which evolves according to x k = Ax k − + Bu k , (10) in the sense of [9, 11] © 2021 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating newcollective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. here k ∈ { , , . . . , K } is the time index (with finite horizon K ), and where both A ∈ R N × N and B ∈ R N × are assumedto be known. Our goal is to determine a two-level input signal u , . . . , u K ∈ { a, b } such that some output (or feature) y k = Cx k ∈ R L (11)(with known C ∈ R L × N ) follows a given trajectory ˘ y , . . . , ˘ y K ∈ R L , i.e., we wish K (cid:88) k =1 (cid:107) y k − ˘ y k (cid:107) (12)to be as small as possible. For ease of exposition, we willassume that the initial state x is known.Note that this offline control problem may be viewed asa single episode of an online control problem, with planninghorizon K . Note also that we are primarily interested in K (cid:29) , which precludes exhaustive tree search algorithms. In order to solve the problem stated in Section 3.1, we turn itinto a statistical estimation problem with an (improper) i.i.d.prior ρ ( u, θ ) (cid:44) K (cid:89) k =1 ρ ( u k , θ k ) , (13)where u (cid:44) ( u , . . . , u K ) , θ (cid:44) ( θ , . . . , θ K ) , and ρ ( u k , θ k ) (cid:44) N (cid:0) u k ; a, σ ,k (cid:1) N (cid:0) u k ; b, σ ,k (cid:1) (14)with θ k = ( σ ,k , σ ,k ) as in (1). Accordingly, we replace (12)by the likelihood function p (˘ y | u ) = K (cid:89) k =1 π ) L/ s L exp (cid:18) −(cid:107) y k − ˘ y k (cid:107) s (cid:19) , (15)where ˘ y = (˘ y , . . . , ˘ y K ) and where s > is a free parameter.The complete statistical model is then given by p (˘ y, u, θ ) (cid:44) p (˘ y | u ) ρ ( u, θ ) (16)together with (10) and (11). Both joint MAP estimation of u and θ (as in Section 2.1) andtype-II estimation of u and θ (as in Section 2.2) can be imple-mented as special cases (with different versions of Step 2) ofthe following algorithm. In (13), ρ is used for two different functions (with different arguments). The algorithm estimates θ and u by alternating the followingtwo steps for i = 1 , , , . . . :1. For fixed θ = θ ( i − , compute the posterior means ˆ u ( i ) k of u k (for k = 1 , . . . , K ) and, if necessary, the poste-rior variances V ( i ) U k of u k , with respect to the probabilitydistribution p ( u | ˘ y, θ ) .2. From these means and variances, determine new NUVparameters θ ( i ) .Note that Step 1 operates with a standard linear Gaussianmodel. In consequence, the required means and variances canbe computed by standard Kalman-type recursions or, equiva-lently, by forward-backward Gaussian message passing, witha complexity that is linear in K .A preferred such algorithm is MBF message passing as in[16, Section V], which amounts to Modified Bryson–Fraziersmoothing [18] augmented with input signal estimation. Thisalgorithm requires no matrix inversion and is numericallyquite stable. θ and u by Joint MAP Estimation In this case, we wish to compute the estimate ˆ u = argmax u max θ p (˘ y | u ) ρ ( u, θ ) . (17)The double maximization (over u and θ ) is naturally imple-mented by alternating maximization, which can be imple-mented by the IKIE algorithm, with Step 2 given by (cid:0) σ ,k (cid:1) ( i ) = (cid:0) ˆ u ( i ) k − a (cid:1) and (cid:0) σ ,k (cid:1) ( i ) = (cid:0) ˆ u ( i ) k − b (cid:1) . (18) In this case, we wish to compute the estimate ˆ θ = argmax θ (cid:90) p (˘ y | u ) ρ ( u, θ ) d u, (19)which can be carried out by expectation maximization [19]with hidden variable(s) u . The update step for θ is θ ( i ) = argmax θ E [log p (˘ y | U ) ρ ( U, θ )] (20) = argmax θ E [log ρ ( U, θ )] , (21)where the expectation is with respect to p ( u | ˘ y, θ ( i − ) . Theupdate (20) turns out to be computable by the IKIE algorithmwith Step 2 given by (cid:0) σ ,k (cid:1) ( i ) = V ( i ) U k + (cid:0) ˆ u ( i ) k − a (cid:1) and (22) (cid:0) σ ,k (cid:1) ( i ) = V ( i ) U k + (cid:0) ˆ u ( i ) k − b (cid:1) . (23) This is obvious for L = 1 . For L > , a little adaptation is required. © 2021 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating newcollective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. .
51 ˘ yy u Fig. 4 . Digital-to-analog conversion as in Section 4.1 withtarget waveform ˘ y (dashed), digital control signal u (bottom),and filter output signal y (solid blue).
1. The algorithms of this section normally converge to alocal (not the global) maximum of (17) or (19).2. The parameter s controls the error (12). If s is chosentoo small, the algorithm may return a nonbinary u .3. Type-II estimation empirically works better than jointMAP estimation (in agreement with [20]). The numer-ical results in Figs. 4 and 5 are obtained with type-IIestimation.
4. EXAMPLES4.1. Digital-to-Analog Conversion
One method for digital-to-analog conversion is to feed acontinuous-time analog linear filter directly with a binaryoutput signal u of a digital processor. This method requiresan algorithm to compute a suitable binary signal u such thatthe analog filter output approximates the desired analog sig-nal ˘ y . A standard approach is to compute u by a delta-sigmamodulator [21], which requires the analog filter to approxi-mate an ideal low-pass filter. By contrast, the method of thispaper works also with simpler (i.e., less expensive) analogfilters.For the following numerical example, the analog filter is a3rd-order low-pass, resulting in the discrete-time state spacemodel A = . − . − . . . − . . . , (24) B = (cid:2) . (cid:3) T , and C = (cid:2) . (cid:3) . Thenumerical results in Fig. 4 are obtained with a = 0 , b = 1 , K = 450 and s = 0 . . .
51 ˘ yy u Fig. 5 . Flappy bird control with check points ˘ y , binary controlsignal u (bottom), and resulting trajectory y (solid blue). The following control problem is a version of the flappy bird computer game [22]. Consider an analog physical systemconsisting of a point mass m moving forward (left to right inFig. 5) with constant horizontal velocity and “falling” verti-cally with constant acceleration g . The { , } -valued controlsignal u affects the system only if u k = 1 , in which case afixed value is added to the vertical impulse. We wish to steerthe point mass such that it passes approximately through asequence of check points, as illustrated in Fig. 5.For this example, we need a slight generalization of (10)–(12) as follows. The state x k ∈ R (comprising the verticalposition and the vertical speed) evolves according to x k = (cid:20) T (cid:21) x k − + (cid:20) /m (cid:21) u k + (cid:20) − T g (cid:21) , (25)and we wish the vertical position y k = (cid:2) (cid:3) x k to minimize K (cid:88) k =1 w k ( y k − ˘ y k ) , (26)where w k = 1 if ˘ y k is a checkpoint and w k = 0 otherwise.The numerical results in Fig. 5 are obtained with m = 0 . , T = 0 . , g = 0 . , a = 0 , b = 1 , K = 250 , and s = 0 . .
5. CONCLUSION
We have proposed a new method for controlling a linear sys-tem with binary input, which can also be used for digital-to-analog conversion. The key idea is a new binary-enforcingprior with a NUV representation, which turns the actual com-putations into iterations of Kalman-type forward-backwardrecursions. The computational complexity of the proposedmethod is linear in the planning horizon, with compares fa-vorably with existing “optimal” methods.The proposed prior and method can be extended both to
M > levels and to sparse level switching in the controlsignal, as will be detailed elsewhere.The suitability of the proposed prior for other applicationsremains to be investigated. This generalization is effortlessly handled by IKIE. © 2021 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating newcollective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. . REFERENCES [1] A. H. Land and A. G. Doig, “An automatic methodof solving discrete programming problems,”
Economet-rica , vol. 28, no. 3, pp. 497–520, 1960.[2] L. A. Wolsey and G. L. Nemhauser,
Integer and Combi-natorial Optimization . John Wiley & Sons, 1999.[3] R. P. Aguilera and D. E. Quevedo, “On stability and per-formance of finite control set MPC for power convert-ers,” in
IEEE Workshop on Predictive Control of Elec-trical Drives and Power Electronics , 2011, pp. 55–62.[4] T. Geyer and D. E. Quevedo, “Multistep finite con-trol set model predictive control for power electronics,”
IEEE Trans. Power Electron. , vol. 29, no. 12, pp. 6836–6846, 2014.[5] T. Dorfling, H. du Toit Mouton, T. Geyer, and P. Kara-manakos, “Long-horizon finite-control-set model pre-dictive control with nonrecursive sphere decoding on anFPGA,”
IEEE Trans. Power Electron. , vol. 35, no. 7, pp.7520–7531, 2020.[6] S. Sparrer and R. F. H. Fischer, “Adapting compressedsensing algorithms to discrete sparse signals,” in , 2014,pp. 1–8.[7] S. Russel and P. Norvig,
Artificial intelligence: A mod-ern approach . Pearson Education Limited, 2013.[8] M. Pincus, “Letter to the editor – a Monte Carlo methodfor the approximate solution of certain types of con-strained optimization problems,”
Operations Research ,vol. 18, no. 6, pp. 1225–1228, 1970.[9] M. E. Tipping, “Sparse Bayesian learning and the rel-evance vector machine,”
Journal of Machine LearningResearch , vol. 1, pp. 211–244, 2001.[10] M. E. Tipping and A. C. Faul, “Fast marginal likelihoodmaximisation for sparse Bayesian models,” in
Proc. ofthe Ninth International Workshop on Artificial Intelli-gence and Statistics , 2003, pp. 3–6.[11] D. P. Wipf and B. D. Rao, “Sparse Bayesian learning forbasis selection,”
IEEE Trans. Signal Process. , vol. 52,no. 8, pp. 2153–2164, 2004. [12] D. P. Wipf and S. S. Nagarajan, “A new view of auto-matic relevance determination,” in
Advances in NeuralInformation Processing Systems , 2008, pp. 1625–1632.[13] H.-A. Loeliger, B. Ma, H. Malmberg, and F. Wadehn,“Factor graphs with NUV priors and iterativelyreweighted descent for sparse least squares and more,”in
Proc. Int. Symp. Turbo Codes & Iterative Inform. Pro-cess. (ISTC) , 2018, pp. 1–5.[14] F. Bach, R. Jenatton, J. Mairal, and G. Obozinski, “Opti-mization with sparsity-inducing penalties,”
Foundationsand Trends in Machine Learning , vol. 4, no. 1, pp. 1–106, 2012.[15] J. Dai, A. Liu, and H. C. So, “Sparse Bayesian learn-ing approach for discrete signal reconstruction,” 2019,unpublished, arXiv:1906.00309.[16] H.-A. Loeliger, L. Bruderer, H. Malmberg, F. Wadehn,and N. Zalmai, “On sparsity by NUV-EM, Gaussianmessage passing, and Kalman smoothing,” in
Informa-tion Theory and Applications Workshop (ITA) , La Jolla,CA, 2016, pp. 1–10.[17] H.-A. Loeliger, “An introduction to factor graphs,”
IEEESignal Process. Mag. , vol. 21, no. 1, pp. 28–41, 2004.[18] G. J. Bierman,
Factorization Methods for Discrete Se-quential Estimation . Academic Press, 1977, vol. 128.[19] P. Stoica and Y. Sel´en, “Cyclic minimizers, majoriza-tion techniques, and the expectation-maximization al-gorithm: a refresher,”
IEEE Signal Proc. Mag. , vol. 21,no. 1, pp. 112–114, 2004.[20] R. Giri and B. Rao, “Type I and type II Bayesian meth-ods for sparse signal recovery using scale mixtures,”
IEEE Trans. on Signal Process. , vol. 64, no. 13, pp.3418–3428, 2016.[21] B. E. Boser and B. A. Wooley, “The design of sigma-delta modulation analog-to-digital converters,”
IEEE J.Solid-State Circuits , vol. 23, no. 6, pp. 1298–1308,1988.[22] Flappy Bird. Accessed 09-October-2020. [Online].Available: https://en.wikipedia.org/wiki/Flappy Bird, vol. 23, no. 6, pp. 1298–1308,1988.[22] Flappy Bird. Accessed 09-October-2020. [Online].Available: https://en.wikipedia.org/wiki/Flappy Bird