Towards Analog Memristive Controllers

1 Towards Analog Memristive Controllers Index Terms—BMI Optimization, Control of Miniaturized Devices, Gain Scheduling, Memristor, Robust and Adaptive Control I. I NTRODUCTION EMRISTOR [1], considered as the fourth basic circuit element, remained dormant for four decades until the accidental discovery of memristance in nanoscopic crossbar arrays by a group of HP researchers [2]. Memristor, an acronym for memory-resistor, has the capability of memorizing its history even after it is powered off. This property makes it a desirable candidate for designing high density non-volatile memory [3]. However, optimal design of such hardware architectures will require accurate knowledge of the nonlinear memristor dynamics. Hence, considerable effort has been channeled to mathematically model memristor dynamics (see [2], [4], [5], [6]). The memorizing ability of memristor has lead researchers to think about its possible use in neuromorphic engineering. Memristors can be used to make dense neural synapses [7] which may find applications in neural networks [8], character recognition [9], emulating evolutionary learning (like that of Amoeba [10]). Other interesting application of memristor may include its use in generating higher frequency harmonics which can be used in nonlinear optics [11] and its use in making programmable analog circuits [12], [13]. Memristor is slowly attracting the attention of the control community. Two broad areas have received attention: 1) Control of Memristive Systems. This include works reported in [14], [15] give detailed insight into modelling and control of memristive systems in Port-Hamiltonian framework while State-of-the-Art work may include [16] which studies global stability of Memristive Neural Networks. 2) Control using Memristive Systems. The very first work in this genre was reported in [17] where the author derived the describing function of a memristor which can be used to study the existence M The authors are with the Department of Electrical Engineering, Indian Institute of Technology, Madras, India. e-mail: {ee13s005, ramkrishna, ee13s007}@ee.iitm.ac.in Voltage V Current I dV / dI = R dt = dΦ / dI = L A / V t= I dΦ w V Charge Q D ? dΦ / dQ = M + Undoped Doped /d dQ / dV = C Abstract—Memristors, initially introduced in the 1970s, have received increased attention upon successful synthesis in 2008. Considerable work has been done on its modeling and applications in specific areas, however, very little is known on the potential of memristors for control applications. Being nanoscopic variable resistors, one can consider its use in making variable gain amplifiers which in turn can implement gain scheduled control algorithms. The main contribution of this paper is the development of a generic memristive analog gain control framework and theoretic foundation of a gain scheduled robustadaptive control strategy which can be implemented using this framework. Analog memristive controllers may find applications in control of large array of miniaturized devices where robust and adaptive control is needed due to parameter uncertainty and ageing issues. dQ arXiv:1405.7677v3 [math.OC] 31 Jul 2014 Gourav Saha, Ramkrishna Pasumarthy, and Prathamesh Khatavkar _ a) Flux Φ + Ron Dw  _ Rof f 1 − Dw b) Figure 1. a) Memristor as a missing link of Electromagnetic Theory. Adapted from [4]. b) Memristor as a series resistor formed by doped and undoped regions of T iO2 . Adapted from [2]. of undesirable limit cycles (i.e. sustained oscillations) in a closed loop system consisting of memristive controller and linear plant. Another line of work may include [18], [19] which uses memristor as a passive element to inject variable damping thus ensuring better transient response. In this paper we lay the groundwork to use memristor as an analog gain control (AGC) element for robust-adaptive control of miniaturized systems. Why "Analog Memristive Controller"?: Several applications needs controlling an array of miniaturized devices. Such devices demands Robust-Adaptive Control due to parameter uncertainty (caused by design errors) and time varying nature (caused by ageing effect). Robust-Adaptive Control algorithms found in literature are so complex that they require micrcontroller for implementation. This poses scalability and integration issues [20] (page 190) because microcontroller in itself is a complicated device. The motive of this work is two folded: 1) Invoke a thought provoking question: “Can modern control laws (like Robust-Adaptive control) be implemented using analog circuits1 ?”. 2) Suggest memristor as an AGC element2 for implementing Robust-Adaptive Control. The paper is organized as follows. We begin our study by gaining basic understanding of memristor in Section II. A generic gain control architecture using memristor is proposed in Section III. Section IV discusses a robust-adaptive control strategy which can be implemented in an analog framework. Section V deals with designing an analog controller for a miniaturized setup by using results from Section III and IV. 1 A microcontroller (or a few microcontrollers) may be used to control an array of miniaturized devices by using Time Multiplexing. In time multiplexing a microcontroller controls each element of the array in a cyclic fashion. In such a scenario a microcontroller will face a huge computational load dependent on the size of the arrray. Upcoming ideas like “event-based control” promises to reduce the computational load by allowing aperodic sampling. The motive of this work is not to challenge an existing idea but to propose an alternative one. 2 CMOS based hybrid circuits like that proposed in [21] can also act like variable gain control element. However memristors are much smaller (found below 10 nm) than such hybrid circuits and hence ensures better scalability. Gain-Span of memristor is also more than CMOS based hybrid circuits. 2 W indow F unction f D !w" 1 0.95 0.8 0.6 0.4 Effective Normalized Channel Length =(0.91−0.08) =0.83 0.2 0 0 wl D 0.08 0.2 0.4 0.6 N ormalized Channel Length 0.8 0.91 w D 1 wh D Figure 2. Plot of jogelkar window function with p = 8 for HP memristor. For
w w the safe zone of operation is D ∈ [0.08 , 0.91]. a tolerance of 5% in f D II. M EMRISTOR P RELIMINARIES Chua [1] postulated the existence of memristor as a missing link in the formulation of electromagnetism. Fig. 1a) gives an intuitive explanation of the same. As evident from Fig. 1a), such an element will link charge Q and flux Φ, i.e. Φ = f (Q). Differentiating this relation using chain rule and applying Lenz’s Law yields V = M (Q) I (1) suggesting that the element will act like a charge controlled df as the variable resistance. variable resistor with M (Q) = dQ This device has non-volatile memory ([1], [4]) and are hence called memristors. Memristive systems [22] are generalization of memristor with state space representation, Ẇ = F (W , I) ; V = R (W , I) I (2) where, W are the internal state variables, I is the input current, V is the output voltage, R (W , I) is the memristance. Memristor reported by HP Labs in [2] is essentially a memristive system with the following state space representation, ẇ = µ RDon f w D
 w I; V = Ron D + Rof f 1 − w D
 I (3) It consists of two layers: a undoped region of oxygen-rich T iO2 and doped region of oxygen-deficient T iO2−x . Doped region has channel length w and low resistance while undoped region has channel length D − w and high resistance. These regions form a series resistance as shown in Fig. 1b. Ron and Rof f are the effective resistance when w = D and w = 0 respectively with Rof f ≫ Ron . µ is the ion mobility. When external bias is applied, the boundary between the two regions drift. This drift is slower near the boundaries, i.e. ẇ → 0 as w → 0 or w → D. This
nature is captured in [4], [5] w . A
plot of Joglekar window using window
function f D p 2w w function f D = 1 − D − 1 is shown in Fig. 2. A close investigation of various proposed models [2], [4], [5],
w [6] reveals two important facts: 1) As shown in Fig. 2, f D is approximately 1, except near the boundaries. 2) Boundary dynamics of memristor is highly non-linear and still a matter of debate. Hence the region w ∈ [wl , wh ], 0 < wl < wh < D,
w ≈ 1 is the safe zone in which memristor where f D dynamics can be approximated as  S  S Q̇M = I ; V = Rof (4) f − α QM I S (RS f −Ron ) S where, αS = ofQ , Rof f = Rof f − (Rof f − Ron ) wDl , S M S Ron = Rof f − (Rof f − Ron ) wDh , DS = wh − wl , QSM = DS D µRon . Superscript “S” means “safe”. In equation (4) we define S QM such that QM = 0 when w = wl . Then QM = Q   M when w = wh . Hence equation (4) is valid when QM ∈ 0 , QSM . From now on, the following conventions will be used: 1) Memristor would mean a HP memristor operating in safe zone. Memristor dynamics is governed by equation (4). 2) The schematic symbol of the memristor shown in Fig. 1a will be used. Conventionally, resistance of memristor decreases as the current enters from the port marked "+". 3) HP memristor parameters: Ron = 100Ω, Rof f = 16kΩ, 2 D = 10nm, µ = 10−14 m sV . We model the memristor using Joglekar Window Function with p = 8. The safe zone is marked by wl = (see
0.08D and wh = 0.91D w S Fig. 2) in which f D ≥ 0.95. This gives: Rof f = Ω S 15kΩ, Ron = 1.5kΩ, QSM = 83µC, αS = 1.6 × 108 C . These parameters will used for design purposes. III. A NALOG G AIN C ONTROL F RAMEWORK In this section we design a AGC circuit whose input-output relation is governed by the following equation, K̇ = αk VC (t) ; Vu (t) = KVe (t) (5) We assume that K > 0. This is basically a variable gain proportional controller with output voltage Vu , error voltage Ve and variable gain K. K is controlled by control voltage VC . αk determines the sensitivity of VC on K. An analog circuit following equation (5) is generic in the sense that it can be used to implement any gain scheduled control algorithm. We assume VC and Ve are band limited, i.e. - the M maximum frequency component of VC and Ve are ωC and M M M ωe respectively. Knowledge of ωC and ωe is assumed. The proposed circuit is shown in Fig. 3. We assume the availability of positive and negative power supply3 VDD and −VDD . All op-amps are powered by VDD and −VDD . Claim 1: The proposed circuit shown in Fig. 3. is an approximate analog realization of equation (5) if: 1) Electrical components in Fig. 3. are ideal. Also for MOSFET’s, threshold 0.  Mvoltage
Vth ≈ M + ωeM , 2ωC 2) ωm = 1000 max ωC 1 3) Rf Cf = 100 ωm ; Re Ce = 2(ω M +ωeM ) C We understand the working of the proposed circuit by studying its four distinct blocks and in the process prove the above claim. The output response of each block for a given input, Ve and VC , is shown in Fig. 4. It should be noted that the tuning rules proposed in the claim is only one such set of values which will make the circuit follow equation (5). Remark 1: Substrate of all NMOS and PMOS are connected to −VDD and VDD respectively. In ON state4 , voltages between −VDD to (VDD − Vth ) will pass through5 NMOS and 3V DD and −VDD are the highest and the lowest potential available in the circuit. From control standpoint, it imposes bounds on control input Vu . 4 With a slight abuse of terminology, a NMOS and a PMOS is said to be in ON state when its gate voltage is VDD and −VDD respectively. 5 A voltage is said to pass through a MOSFET if the exact voltage applied at its source(drain) terminal appears at its drain(source) terminal. 3 VL1 Rs Integrator _ VIG T1 Cf Rf Vm + + RI − M Rf D1 + − RC + D2 TC1 Ce Re _ + High Pass Filter Envelope Detector VED Vu − Vf sin (ωm t) Memristor Gain Block − CH VC VCm + − D4 TC4 O2 O1 Cs + D3 T4 T2 + T3 − TC3 TC2 VL2 Ve B. High Pass Filter (HPF) Charge Saturator Block Voltage Follower Figure 3. Memristive AGC Circuit. Ve and VC are the inputs. Vu , VL1 and VL2 are the outputs. VL1 and VL2 are zone indicating voltages. voltages between − (VDD − Vth ) to VDD will pass through PMOS. If Vth ≈ 0, any voltage between −VDD to VDD will pass through NMOS and PMOS when they are in ON state. A. Memristor Gain Block The key idea of this block has been adapted from [23]. VCm and Ve are inputs to this block while Vm is the output. For now we assume VCm = VC (details discussed in III-D). Current Im through memristor is Im (t) = Ve (t) sin (ωm t) VC (t) + RI R {z } | {zC } | Inoise Icntrl From equation (4), resistance of the memristor is given by S S M (QM ) = Rof f − α QM . Differentiating this relation we S get, Ṁ = −α Im (t). Hence voltage Vm is given by Ṁ = −αS Im (t) ; Vm (t) = −Im (t) M (6)  M
M Note that ωm = 1000 max ωC + ωeM , 2ωC > ωeM . In such a case the minimum component frequency of Inoise = Ve (t) sin(ωm t) is ωm − ωeM . Now RI  M
M ωm − ωeM = 1000 max ωC + ωeM , 2ωC − ωeM
M ≥ 1000 ωC + ωeM − ωeM = M M 1000ωC + 999ωeM ≫ ωC This implies that the lowest component frequency of Inoise is much greater than the highest component frequency of Icntrl . Also note that Ṁ = −αS Im (t) in an integrator with input Im (t) and output M (t). As integrator is a low pass filter, the effect of high frequency component Inoise on M is negligible compared to Icntrl . Hence equation (6) can be modified as αS VC (t) ; Vm (t) = Vm1 + Vm2 where, (7) Ṁ ≈ − RC Vm1 Ve (t) − MR I Vm2 VC (t) . − MR C Vm1 is M . =R I sin (ωm t) and Note = = the modulated form of the desired output with gain K Remark 2: M (t) is the variable gain. According to equation (5), Ve should not effect M (t). Without modulating Ve , the effect of Ve (t) on M (t) would not have been negligible. The role of HPF is to offer negligible attenuation to Vm1 and high attenuation to Vm2 thereby ensuring that the Envelope Detector can recover the desired output. Note that M (t) is basically the integral of VC (t). Since integration is a linear operation it does not do frequency translation. Hence the component frequencies of M (t) and VC (t) are same. Let ω1m and ω2M denote the minimum component frequency of Vm1 and maximum component frequency of Vm2 respectively. Now ω1m = = = ω2M ωm − Maximum Frequency Component of M Ve ωm − Maximum Frequency Component of VC Ve
M ωm − ωC + ωeM (8) = Maximum Frequency Component of M VC = Maximum Frequency Component of (VC ) = M 2ωC 2 (9) The  attenuation offered  by the HPF is given by −2 dB. We want to study the atten−10 log 1 + (ωRf Cf ) uation characteristics at two frequencies: 1) At ω = ω1m : At this frequency ω1m Rf Cf

  M M 100 ωm − ωC + ωeM ωC + ωeM = 100 1 − = ωm ωm   1 ≈ 100 (10) ≥ 100 1 − 1000 Inequality because ωm =  M(10) M
is possibe M m implying ωMω+ω 1000 max ωC + ωe , 2ωC M ≥ 1000. e C For ω1m Rf Cf ≥ 100 attenuation is approximately 0 dB. As the HPF offers almost no attenuation to the minimum component frequency of Vm1 , it will offer no attenuation to the higher component frequency of Vm1 as well. Hence Vm1 suffers negligible attenuation. 2) At ω = ω2M : At this frequency M 100 2ωC ≤ = 0.1 (11) ωm 1000 Inequality because ωm =  M(11) M
is possibe ωm M implying 2ω ≥ 1000. 1000 max ωC + ωe , 2ωC M C For ω2M Rf Cf ≤ 0.1 attenuation is more than −20 dB. As the HPF offers high attenuation to the maximum component frequency of Vm2 , it will offer higher attenuation to the lower component frequency of Vm2 . Hence Vm2 gets highly attenuated. As Vm1 undergoes almost no attenuation (≈ 0dB), the output of HPF is Vf = −Vm1 . The minus sign before Vm1 is justified as the HPF is in inverting mode. ω2M Rf Cf = 100 C. Envelope Detector Ve (t) sin (ωm t). Similar The input to this block is Vf = M R I to amplitude modulation (AM), here sin (ωm t) is the carrier Ve (t) and M R is the signal to be recovered. We use a polarity I Ve (t) sensitive envelope detector as M R can be both positive or I negative. The key idea used here is that the polarity of Ve and C e modulated, i.e. MRVI e , may  contain frequencies anywhere in M M the range 0, ωC + ωe . For a given Re Ce a signal with a lower frequency will suffer higher ripple. Hence to constrain the ripple for any frequency in the specified range one must constrain the ripple for 0 frequency (DC voltage). For DC 2π×100 . With the voltage ripple factor is approximately √3ω m R e Ce choice of Re Ce made in Claim 1 and 100 as multipying factor, ripple factor is as large as 7.2%. Therefore, we choose multiplying factor of 1000 which gives a ripple factor of Ve (t) 0.72%. The output of the envelope detector is M R and I the final output of the circuit is Ṁ ≈ − where, M = αS VC (t) ; Vu (t) = MVe RI R C (12) M RI . Comparing equations (5) and   (12) we see S S Rof Ron f α that αk = − RI RC and K = M where M ∈ RI , RI . S RI is tuned to get the desired range of gain while RC is a free parameter which can be tuned according to the needs. D. Charge Saturator Block This block limits the memristor to work in its safe zone hence ensuring validity of equation (4). In safe zone the Vm (in V ) 0.5 0.25 0 −0.25 −0.5 5 2.5 0 −2.5 0 200 400 600 800 a) t (in ms) −5 1000 5 0 −5 −10 0 200 400 600 800 200 400 600 800 1000 0 200 400 600 800 1000 Vu (in V ) 0 −1 200 400 600 800 c) t (in ms) 0 200 400 600 800 1000 0 −1 e) t (in ms) 2 1 0 d) t (in ms) 1 −2 1000 b) t (in ms) 2 −2 0 2 Vf (in V ) VC (in V ) 10 True Output (in V) M > 0. Hence we detect the positive Vu is same since K = R I peaks of Vf when Ve is positive by keeping T1 ON and TC1 OF F . When Ve is negative, negative peaks of Vf are detected by keeping TC1 ON and T1 OF F . Remaining working of the envelope detector is similar to a conventional Diode-based Envelope Detector and can be found in [24]. Effective envelope detection using diode based envelope detector requires: M 1) ωm ≫ ωC + ωeM , i.e. the frequency of the carrier should be much greater than the maximum component frequency of the signal. Here sin (ωm t) is the carrier whose frequency is ωm while MRVI e is the signal whose M maximum component frequency is ωC + ωeM (refer equation (8)). 2) ω1m ≪ Re Ce , i.e. the time constant of the envelope detector is much larger than the time period of the modulating signal. This ensures that the capacitor Ce (refer Fig. 3) discharges slowly between peaks thereby reducing the ripples. 1 3) Re Ce < ωM +ω M , i.e. the time constant of the envelope e C detector should be less than the time period corresponding to the maximum component frequency of the signal getting modulated. This is necessary so that the output of the envelope detector can effectively track the envelope of the modulated signal. The proposed tuning rule in Claim 1 satisfies these conditions. In general, for amplitude modulation (AM) the choice of modulating frequency is 100 times the maximum component frequency of the signal getting modulated. Unlike AM, our multiplying factor is 1000 instead of 100. The reason for choosing this should be clearly understood. In a conventional diode based peak detector the ripple in output voltage can be decreased either by increasing the modulating frequency ωm or by increasing the time constant Re Ce . But Re Ce is 1 upper bounded by ωM +ω M . In our case the signal to be Ve (in V ) 4 1000 1 0 −1 −2 f) t (in ms) Figure 4. Output of various stages of the circuit shown in Fig. 3, i.e. Vm , Vf and Vu , corresponding to inputs: Ve and VC . Parameters of simulations are: HP Memristor in safe zone, ωm = 628 × 103 rads−1 , RC = 100 kΩ, RI = 1 kΩ, Rf Cf = 0.159 ms, Re Ce = 0.796 ms, VDD = 5 V , Rs Cs = 0.826 s. Vu is obtained by simulating the circuit show in Fig. 3. The true output is obtained by numerically solving equation (5). In both cases we use Ve and VC shown in Fig. 4a and Fig. 4b respectively as inputs. following equations are valid6 , dVIG dVIG RC V m dQM Vm =− (13) =− C ; = C ⇒ dt Rs C s dt RC dQM Rs C s  that in the safe zone w ∈ [wl , wh ] and QM ∈  Recall 0 , QSM . We assume that integrator voltage VIG = 0 when QM = 0 (or w = wl ). Integrating equation (13) under this assumption yields RC QM (14) VIG = − Rs C s In equation (14),hif QM = QSM i (or w = wh ), VIG = − R QS − RCs CM s RC QS M R s Cs . Hence, VIG ∈ , 0 when memristor is in its safe zone. In Fig. 3 comparator O1 and O2 are used to compare VIG to know if the memristor is in safe zone. Note that: 1) Reference voltage of comparator O1 and O2 is GN D and voltage R QS VH (across capacitor CH ) respectively. VH is set to − RCs CM s by Synchronization Block (refer Section III-E). 2) In Fig. 3 any MOSFET transistor couple, Ti and TCi , are in complementary state, i.e. if Ti is ON , TCi will be OF F and vice-versa. 3) Comparator output VL1 and VL2 gives knowledge about the state of the memristor. Also VL1 , VL2 ∈ {−VDD , VDD }. Now depending on VL1 and VL2 , three different cases may arise: Case 1 (VH < VIG < 0 ⇒ VL1 = VDD , VL2 = VDD ) This implies that the memristor is in its safe zone, i.e. 0 < QM < QSM . QM can either increase or decrease. Hence both T2 and T3 are ON allowing both positive and negative voltage to pass through, i.e. VCm = VC . Case 2 (VIG < VH ⇒ VL1 = VDD , VL2 = −VDD ) This happens when QM ≥ QSM . Since QM can only decrease, T2 is kept OF F but T3 is ON . Two cases are possible: if 6 Actually, dQM dt Vm V sin((ω t)) = RC + e R m . But as proved in III-A, the effect C I of high frequency term, Ve sin (ωm t) on QM is negligible. 5 A C M _ TC6 T6 RC VL3_ D VL3 15 −1 10 V 1 5 V −2 _ O3 S S3 Ron VDD 20 −1.5 V1 + 0 −0.5 V 1 , V 2 (in V ) VL3 T5 + TC5 M (in kΩ) B VL3 S S1 Rof f VL3 + 2 0 0.5 1 1.5 0 2 b) t (in s) 12 0 0.5 1 1.5 2 c) t (in s) 0 V2 + VL3 S4 S2 Cs Rs − Integrator VIG V I G (in V ) _ V L3 (in V ) −2 RC 0 Rapid Switching −6 −8 −10 CH −12 + −4 0 0.5 1 1.5 2 −12 0 0.5 a) 1 1.5 2 e) t (in s) d) t (in s) Figure 5. a) Schematic of Synchronization Block. b), c), d), e) Graphs showing operations of Synchronization Block. In these graphs preset mode operates for the first 1s and online calibration mode operates for the next 1s. VC > 0, T4 and TC2 will be ON making VCm = 0 or VC < 0 making it pass through T3 and thereby setting VCm = VC . Case 3 (VIG > 0 ⇒ VL1 = −VDD , VL2 = VDD ) This happens when QM ≤ 0. Since QM can only increase, T3 is kept OF F but T2 is ON . Two cases are possible: if VC < 0, TC4 and TC3 will be ON making VCm = 0 or VC > 0 making it pass through T2 and thereby setting VCm = VC . E. Synchronization Block Operation of Charge Saturator Block assumes that: 1) VIG = 0 when QM = 0 (or w = wl ). 2) Voltage VH across R QS capacitor CH equals − RCs CM . This block ensures these two s conditions and thereby guarantees that the memristor and the integrator in Fig. 3 are in “synchronization”. Synchronization Block is shown in Fig. 5a). In Fig. 5a, the op-amp with the memristor, the integrator and capacitor CH are indeed part of the circuit shown in Fig. 3. Such a change in circuit connection is possible using appropriate switching circuitry. This block operates in two modes: Preset Mode: We first ensure that VIG = 0, when w = wl . In this mode switch S1 and S2 are ON and switch S3 and S4 are OF F . When S2 is closed the residual charge in capacitor CH will get neutralized thus ensuring VIG = 0. Next we make w = RS VDD wl . Note that V1 = − MRVCDD and V2 = − ofRf C . If M < S Rof f then V1 > V2 ⇒ VL3 = −VDD . Hence the path ADBC of wheatstone bridge arrangement will be active making the ve ve current flow from (−) to (+) terminal of the memristor. S S This will increase M till M = Rof f . If M > Rof f , the path S ABDC is active making M decrease till M = Rof f. 7 Online Calibration : Immediately after preset mode is complete, S1 and S2 are switched OF F and S3 and S4 are RS V switched ON . Now, V1 = − MRVCDD and V2 = − onRCDD . S In this step VL3 = VDD as M > Ron always. Path ABDC S will be active driving M to Ron . As VL3 is given as an R QS we tune Rs and Cs s.t. − RC CM = −VDD , VIG ∈ [−VDD , 0] when s s memristor is in safe zone. Then we can directly use power supply −VDD as the reference voltage for comparator O2 . This will eliminate the need 7 If R QS of capacitor CH and “Online Calibration”. However if − RC CM 6= −VDD s s (due to tuning error), this approach may drive the memristor to non-safe zone. input to the integrator, capacitor CH will also get charged. Note that in this step memristor will work in safe zone. Also VIG = 0 when QM = 0 (ensured by preset mode). Hence relation between VIG and QM will be governed by equation S (14). When M gets equalized to Ron , QM = QSM , thereby RC QS . making VH = VIG = − Rs CM s Each of the modes operate for a predefined time. The resistors and hence the voltages V1 and V2 may get equalized before the predefined time after which VL3 will switch rapidly. Such rapid switching can be prevented by replacing the comparator O3 by a cascaded arrangement of a differential amplifier followed by a hysteresis block. Various graphs corresponding to synchronization process are shown in Fig. 5 b), c), d), e). Memristor Gain Block and the Integrator of Charge Saturator Block should be periodically synchronized to account for circuit non-idealities. One such non-ideality can be caused if the capacitor CH is leaky causing the voltage VH to drift with time. Remark 3: In the discussion of the Synchronization Block S S we have slightly misused the symbols Ron and Rof f . ResisS S tance of Ron and Rof f shown in Fig. 5a should be close S S to the actual Ron and Rof f (as mentioned in Section II) respectively. It is not necessary that they should be exactly equal. However, there resistance must lie within the safe zone. This alleviates analog implementation by eliminating the need of precision resistors. It should be noted that the maximum and the minimum resistance of the memristor in safe zone S S is governed by the resistances Ron and Rof f used in the S S Synchronization Block not the actual Ron and Rof f. To conclude, in this section we designed an Analog Gain Control framework using Memristor. Schematic of Memristive AGC is shown in Fig. 2 whose circuit parameters can be tuned using Claim 1. Memristive AGC’s designed in this work is “generic” in the sense that it can be used to implement several Gain-Scheduled control algorithms. IV. C ONTROL S TRATEGY As mentioned in the introduction, we are interested in designing Robust Adaptive Control Algorithms to tackle issues like parameter uncertainty (caused by design errors) and time varying nature (caused by ageing effect and atmospheric variation). ‘Simplicity’ is the key aspect of any control algorithm to be implementable in analog framework as we do not have the flexibility of ‘coding’. Robust Adaptive Control Algorithms found in control literature cannot be implemented in an analog framework due to their complexity. Here we propose a simple gain-scheduled robust adaptive control algorithm which can be easily implemented using Memristive AGC discussed in Section III. We prove the stability of the proposed algorithm using Lyapunov-Like method in Section IV-A. Notations: The notations used in this paper are quite standard. R+ and Rn denotes the set of positive real numbers and the n-dimensional real space respectively. I represents identity matrix of appropriate dimension. |·| is the absolute value operator. ∅ represents a null set. The bold face symbols S and S+ represents the set of all symmetric matrices and positive definite symmetric matrices respectively. inf (·) (sup (·))  Rest to Scan Mode Ė  h  E Scan to Rest Mode  6 Legend RGS Gain kM−km T Stabilizing Gain Set Searching t = t1 −γα 0 Ė E kM km Searching t = t2 a) kM km Case 1 Searching Searching t = t3 t = to Searching t = to + T Found t = to + 2T kM km Searching t = t4 km kM Searching t = t5 Searching kM Case 2 km Found t = to Searching km t = to + T kM c) t = t5 + δt Searching t = t5 + δt Found t = to + 2T Reflection at kM −α km Jump NOT Allowed kM b) Figure 6. a) Hysteresis Function as mentioned in equation (16). T is the scan time. α ∈ R+ − {0} and γ ∈ [0, 1) are tuning constants. b) First six images shows the scan/rest mode of RGS. The last three images depicts the concept of drifting. c) Figure showing the worst case scan time. represents the infimum (supremum) of a set. For a matrix A, λm (A) and λM (A) denotes the minimum and maximum eigenvalue of A respectively. A  B implies B −A is positive semi-definite while A ≺ B implies B − A is positive definite. The euclidean norm of a vector and the induced spectral norm of a matrix is denoted · . The operator × when applied on sets implies the cartesian product of the sets. Conv {A} implies the convex hull of set A. Analysis presented futher in this paper uses ideas from real analysis, linear algebra and convex analysis. For completeness, these ideas are briefly reviewed in Appendix C. A. Reflective Gain Space Search (RGS) Consider a SISO variant (LTV) system  uncertain linear time T with states x = x1 x2 · · · xN ∈ RN , input u ∈ R and output y ∈ R described by the following state space equation ẋ = A (t) x + B (t) u ; y = Cx = x1 where, (15)     0 1 ··· 0 0    ..  .. .. .. ..     . . . . A (t) =   B (t) =  .   0   0  0 ··· 1 a1 (t) a2 (t) · · · aN (t)  b (t) The output matrix C = 1 · · · 0 0 . Assumptions: 1) The order N of the system is known. 2) The variation of ai (t) , ∀i = 1, 2, . . . , N and b (t) due to parameter uncertainty and time variation is bounded, i.e. (A (t) , B (t)) belongs to a bounded set L. L is assumed to be a connected set. A (t) and B (t) are not known but knowledge of L is assumed. This assumption basically means that we don’t have knowledge about ai (t) and b (t) however we know the range in which they lie. 3) Ȧ ≤ δA and Ḃ ≤ δB . Knowledge of δA and δB is assumed. This assumption basically means that A (t) and B (t) are continous in time. 4) Let b (t) ∈ [bl , bu ] s.t. bl ≤ bu , and either bl > 0 or bu < 0. This choice of bl and bu is explained in Section IV-B. Example to clarify the concept of L, δA and δB is dealt later in Section IV-B however we give the following example to better explain Assumption 4. Example 1: Consider two cases: 1) b (t) ∈ [−2, 3] 2) b (t) ∈ [0.5, 3]. According to Assumption 4, Case 2 is possible while Case 1 is not. This example clearly illustrates that according to Assumption 4 the sign of b (t) does not change with time and is known with absolute certainty. Note that sign of b (t) decides the sign of static gain8 of the system. So Assumption 4 in certain sense means that the sign of static gain does not change with time and is known with absolute certainty. For all practical scenario such an assumption is valid. Notice that all assumptions are mild from practical viewpoint. Our aim is to regulate the output around the operating point y ∗ = x∗1 = r. Conventional set-point control consist of a bias term ub (t) plus the regulation control input ur (t), i.e. u (t) = ub (t) + ur (t). Here we assume that ub (t) is  T designed s.t. x∗ = r 0 . . . 0 is the equilibrium point and concentrate on the synthesis of ur (t). For simplicity  T we consider r = 0, i.e. x∗ = 0 0 . . . 0 is the equilibrium point. As we are dealing with a linear system the same analysis is valid for r 6= 0. The controller structure is, ur (t) = −K (t) y (t), where K (t) ∈ [km , kM ] is the variable gain s.t. 0 < km < kM . Let E = xT P x be the Lyapunov Candidate Function9 with P ∈ S+ . Then RGS is as simple Gain-Scheduled control strategy given by the following equation ! Ė ; ur = −Ky (16) K̇ = sgn · h E   is a hysteresis function where, sgn ∈ {−1, 1} and h Ė E shown in Fig. 6a. Working of RGS can be explained as: 1) RGS finds the stabilizing gain10 , i.e. the gain which renders Ė < 0 (Ė < −αE in a more strict sense), by reflecting back and forth between [km , kM ]. RGS is said to be in "Scan Mode" when it is scanning for the stabilizing gain. It goes to "Rest Mode" when stabilizing gain is found. RGS Scan Cycle is clearly depicted in the first six images of Fig. 5b). 2) RGS uses Ė as stability indicator. Ė is found by differentiating E which in turn is calculated using E = T x (t) P x (t). To get the states x (t) we differentiate the output y (t) N − 1 times. 3) Scan Mode is triggered when Ė > −γαE (refer Fig. 6a). Scan Mode is associated with a scan time of T , i.e.  time taken to scan from km to kM . Hence, = kM T−km . The value of sgn is 1 when gain h Ė E space is searched from km to kM and −1 otherwise. Scan mode operates till Ė > −αE. 8 Unlike LTI system, static gain for LTV system may not be well defined. b(t) Here we define static gain of LTV system (15) as a (t) . 1 9 Use of time invariant Lyapunov Function to analyse stability of LTV systems has been used in [25] (Chap. 5, 7) and [26] (Chap. 3). 10 The term ”stabilizing gain” has been slightly misused. Stability of LTV system cannot be assured even if the closed loop system matrix, A (t) − B (t) K (t) C, has negative real eigen part for all t > 0 ([27], [28]). 7 4) Rest is triggered when Ė < −αE. In this mode   Mode Ė = 0, i.e. the stabilizing gain is held constant. h E Rest mode operates till Ė < −γαE. 5) In the process of finding stabilizing gain, LTV system may expand (E increases) in Scan Mode and will contract (E decreases) in Rest Mode. RGS ensures that even in the worst case, contraction is always dominant over expansion, guaranteeing stability. Stabilizing Gain Set11 KP,s (t) for a given P ≻ 0 and s ∈ R+ at time t is defined as n T KP,s (t) = K ∈ R : AC (t) P + P AC (t)  −sI, o AC (t) = A (t) − B (t) KC, km ≤ K ≤ kM Lemma 1: T If KP,s (t) 6= ∅, ∀t ≥ 0 then under Assumption 3, KP,s (t) KP,s (t + δt) 6= ∅ if δt → 0. Proof: This lemma basically means that stabilizing gain set will drift. If the stabilizing gain set is drifting then there will be an intersection between stabilizing gain sets at time t1 and t2 if t2 − t1 is small. Drifting is obvious under Assumption 3 however a formal proof is given in Appendix A. Concept of “drifting” is clearly depicted in the last three images of Fig. 5b. Notice that at t = t5 the stabilizing gain set is almost found. Two possible cases are shown in Fig. 5b for t = t5 + δt where δt → 0. In the first case the stabilizing gain set is drifting while in the other case the stabilizing gain set jumps. In the first case the stabilizing gain is found and RGS goes to Rest Mode. In the second case RGS will just miss the stabilizing gain set. Hence if the stabilizing set keeps jumping, the scan mode may never end. Hence if the stabilizing set keeps jumping, the scan mode may never end. As mentioned in Lemma 1, stabilizing set KP,s (t) never jumps if Assumption 3 is valid. Therefore Lemma 1 is important to guarantee an upper bound on the time period of Scan Mode from which we get the following Lemma. Lemma 2: If KP,s (t) 6= ∅, ∀t ≥ 0 then under Assumption 3, the maximum time period of Scan Mode is 2T . Proof: It is a direct consequence of Lemma 1. The worst case time period of 2T happens only in two cases. Both the cases are shown in Fig. 6c. Theorem 1: LTV system characterized by equation (15) and set L is stable under RGS Control Law proposed in equation (16) if it satisfies the following criterion 1) [C1] If corresponding to all (A, B) ∈ L there exist atleast one gain KAB ∈ [km , kM ] s.t. Proof: Consider the Lyapunov candidate E = xT P x which when differentiated along systems trajectory yields Ė = xT Q (t) x ≤ λM (Q (t)) x T for some s ∈ R+ . 2) [C2] Scan time T < λm (P )α 2 (1−γ ) 2 8δ max(λL M ,0) λL M 11 The stabilizing gain set KP,s (to ) at time t = to , is just a “LyapunovWay” of describing the set of gains which will stabilize the corresponding LTI system (A (to ) , B (to ) , C) at time t = to . ≤ Ė 2 (18) λL M x 2 ≤ λL M E λm (P ) (19) We want to find the maximum possible increase in E in the Scan Mode. At this point it is important to understand that [C1] is another way of saying that KP,s (t) 6= ∅, ∀t ≥ 0. Hence if [C1] is true, then Lemma 2 can be used to guarantee that the maximum period of scan mode is 2T . Let E = Es and t = ts at the end of scan mode. Maximum expansion in E is obtained by integrating inequality (19) from t = t∗ to t = t∗ + 2T ES ≤ βs E ∗ where, (20)  L 2λ T βs = exp λmM(P ) , the worst case expansion factor. The end of scan mode implies that a stabilizing gain is found, i.e. at t = ts , Ė (ts ) ≤ −αE (ts ) (refer Fig. 6a), and the rest mode starts. Here it is important to note the relation between α and s (refer [C1]). [C1] assures that at t = ts a gain K can be found s.t.  Ė (ts ) ≤ λM (Q (ts )) kxk 2 = ≤ 2 −s kxk s E (ts ) (21) − λM (P ) Comparing inequality (21) with Fig. 6a we get α = λMs(P ) . Let τ be the duration of rest mode. We want to find the minimum possible decrease in E in rest mode before it again goes back to scan mode. Again, Ė = xT Q (t) x = xT [Q (ts ) + ∆ (t)] x = Ė (ts ) + xT ∆ (t) x where, (22) ∆ (t) = Q (t) − Q (ts ). As K (t) = K (ts ) ; ts ≤ t ≤ ts + τ (K̇ = 0 in rest mode), ∆ (t) can be expanded as . Here γ ∈ [0, 1), δ = δA + δB kM and = max λM (Q) where the set L Q∈SQ n T L SQ = Q : Q = (A − BKC) P + P (A − BKC) o (A, B) ∈ L, K ∈ [km , kM ] ≤ λL M x where, Q (t) = AC (t) P + P AC (t) and AC (t) = A (t) − B (t) K (t) C. Definitely, Q (t) ∈ S as it is a sum of T the matrix P AC (t) and its transpose AC (t) P . Note that, L λM (Q (t)) ≤ λM , ∀ (A, B) ∈ L and ∀K ∈ [km , kM ]. When λL M < 0, stability is obvious as according to inequality (18) Ė < 0, ∀ (A, B) ∈ L and ∀K ∈ [km , kM ]. As Ė < 0 for any ∀K ∈ [km , kM ], stability if guaranteed even if scanning is infinitely slow, i.e. scan time T → ∞. This is in accordance with [C2]. ∗ We now consider the case when λL M > 0. Say at time t = t , ∗ E = E and the system goes to scan mode to search for a stabilizing gain. From inequality (18) we have, T (A − BKAB C) P + P (A − BKAB C)  −sI (17) 2 ∆ (t) = T [∆A (t) − ∆B (t) K (ts ) C] P +P [∆A (t) − ∆B (t) K (ts ) C] (23) where, ∆A (t) = A (t) − A (ts ) and ∆B (t) = B (t) − B (ts ). Taking norm on both side of equation (23) yields ∆ (t) ≤ ≤ ≤ 2 P [∆A (t) − ∆B (t) K (ts ) C] ∆A (t) + ∆B (t) K (ts ) 2 P 2λM (P ) δ∆t C
(24) 8 where, δ = δA + δB kM and ∆t = t − ts . Note that C = 1. Also kP k = λM (P ) for all P ∈ S+ . Getting back to equation (22) we have Ė = ≤ ≤ ≤ Ė (ts ) + xT ∆ (t) x ≤ −s x −s x 2 + ∆ (t) x 2 + ∆ (t) x 2 Er βr = exp − α 2 (1−γ 2 4δ α (1 − γ) 2δ (26) where, (28)   2 α (1−γ 2 ) 2λL T + λmM(P ) . If β ∈ (0, 1), there β = βs βr = exp − 4δ will be an overall decrease in E in a rest-scan cycle. β ∈ (0, 1) can be assured if
λm (P ) α2 1 − γ 2 (29) T < 8δλL M Now we want to prove that inequality (29) ensures that E (t) → 0 (and hence12 x (t) → 0) as t → ∞.Before proceeding forward we note two things: 1) Theoretically . predictable duration of a rest-scan cycle is Trs = 2T + α(1−γ) 2δ 2) Trs is a conservative estimate of the duration of a rest-scan cycle. Hence it may as well happen that the rest mode lasts for a longer time leading to unpredictable contraction. This is clearly shown in Fig. 7. Now we want to put an upper bound on E (t). Say we want to upper bound E (t) in any one blue dots shown in Fig. 7, i.e. in the green zones. This can be done as follows (30) where, E = Eo at t = 0, η (t) is the number of complete rest-scan cycle13 before time t. In inequality (30), β η(t) is the predictable contraction factor contributed by the red zones in Fig. 7 and the last term is the unpredictable contraction factor contributed by the green zones in Fig. 7. Note that in green zones all we can assure is that, Ė ≤ −γαE, which when integrated yields the last term of inequality (30). Now we want to upper bound E (t) in one of the red dots shown in Fig. 7 which can be done as q E where E = xT P x ≥ λm (P ) kxk2 it implies kxk ≤ λm (P ) λm (P ) > 0 as P ≻ 0. Hence if E → 0 then kxk → 0 implying x → 0. 13 Without additional knowledge of system dynamics, η (t) is not predictable. 12 As tp2 t2 TIME tp3 t3 tp4 0 −γ α −α TIME Figure 7. Plot showing the actual and predictable duration of
rest-scan cycle. Red zones shows the predictable duration, i.e. tpi − ti−1 = Trs , while the actual duration is ti − ti−1 . Also note that EiR ≤ βEiS , the theoretically predictable contraction in a rest-scan cycle. The green zones shows the unpredictable rest modes.   + E (t) ≤ Eo βs β η(t) exp −γα (t − Trs − η (t) Trs ) (31) + βE ∗ E (t) ≤ Eo β η(t) exp (−γα (t − η (t) Trs )) t1 Ė E 0 ≤ βr Es where, (27)  ) , the worst case contraction factor. ≤ tp1 (25) Using inequality (20) and (27) we get Er E2S E2R 0 . Let E = Er at the end of rest mode. Hence, τmin = α(1−γ) 2δ We integrate inequality (25) from t = ts to t = ts + α(1−γ) 2δ to find the minimum possible contraction in rest mode,  E 2 To find the minimum possible value of τ we can substitute Ė = −γαE in inequality (25). This gives, ≥ Unpredictable Rest Modes E1R 2 − [s − 2λM (P ) δ∆t] x − (α − 2δ∆t) E τ Predictable Rest−Scan Cycle Duration E1S In inequality (31), the operator (a) = max (0, a). Inequality (31) resembles (30) except that in this case the system may be in scan mode leading to the extra expansion factor βs in the expression. There is an extra −Trs in the last term to deduct the predictable time of the current rest-scan cycle. Among inequality (30) and (31) , (31) is definitely the worst case upper bound on E (t) due to the presence of two additional terms, i.e. −Trs inside exp (−γα (·)) and βs ≥ 1. From inequality (31) it is clear that if β ∈ (0, 1) then, E (t) → 0 as t → ∞. This completes the proof of Theorem 1. Remark 4: From inequality (31) convergence is better if β is small. According to inequality (28), β decreases as α (or s) increases and T decreases. The effect of γ on β is more involved. In inequality (31), if γ increases then β increases (refer inequality (28)). Thus predictable contraction decreases (due to large β) but unpredictable contraction increases (due to large γ). Hence γ can be neither too high nor too low. Remark 5: Note that RGS strategy doesn’t impose any theoretical limit on the rate of variation, δA and δB , of LTV system. This is perhaps one of the novelty of RGS. Remark 6 (RGS for LTI Systems): For uncertain LTI systems, Theorem 1 will have [C1] without any change. However, [C2] will no longer impose an upper bound on T but will just demand a finite non-zero value. This will ensure that a stabilizing gain is found within a finite time period of 2T . B. Bilinear Matrix Inequality(BMI) Optimization Problem Foundation of RGS is based on the validity of [C1] for a given uncertain LTV/LTI system. We pose a BMI optimization problem to check the validity of [C1] and in the process find the value of P and s needed for implementing RGS. We start our discussion by formally defining the set L. Let A (t) and B (t) be functions of p independent physical T parameters Θ = [Θ1 , Θ2 , . . . , Θp ] , i.e. (A (t) , B (t)) = F (Θ (t)). Θ (t) is time varying and at every time instant is 9 also associated with a uncertainty set because of parameter uncertainty. We assume that every  physical parameter Θi is H L Θ , Θ bounded, i.e. Θ (t) ∈ Θ i i i . Then   L H  (t) ∈ S, where  L H H S = ΘL 1 , Θ1 × Θ2 , Θ2 ×. . .× Θp , Θp a p−dimensional hypercube. We assume no knowledge of the time variation of Θ, i.e. Θ (t), but we assume the knowledge of S. Then L is the image of S under the transformation F , i.e. F : S → L where S ⊂ Rp and L ⊂ RN ×N × RN ×1 . Remark 7: The system described by equation (15) can be represented   using the compact notation a1 a2 · · · aN | b . Hence one can assume that L ⊂ RN +1 rather than L ⊂ RN ×N × RN ×1 . This reduces notational complexity by making the elements of L a vector rather than ordered pair of two matrix. At this point we would be interested in formulating [C1] as an optimization problem. With a slight misuse of variable s we can state the following problem. Problem 1: [C1] holds if and only if the optimal value of the problem minimize: s subject to: T (A − BKAB C) P +P (A − BKAB C)  sI ∀ (A, B) ∈ L P ≻0 with the design variables s ∈ R, P ∈ RN ×N and KAB ∈ [km , kM ] is negative. Note the use of KAB instead of K in Problem 1. It is to signify that we do not need a common gain K for all (A, B) ∈ L. Perhaps we can have seperate gains KAB for every (A, B) ∈ L satisfying Problem 1 and the RGS strategy described in Section IV-A will search for it. However the optimization problem described in Problem 1 is semi-infinite14 and is hence not computationally tractable. We will now pose Problem 1 as a finite optimization problem. We can always bound the compact set L by a convex polytope. Define a polytope15 P = Conv {V} where V = {(Ai , Bi ) : i = 1, 2, . . . , m}, the m vertices of the convex polytope s.t., if (A, B) ∈ L then (A, B) ∈ P. Then L ⊆ P. We now give the following example to illustrate concepts related to S, L and P (discussed above) and also discuss how to calculate δA and δB (discussed in Assumption 3). Example 2: Consider the following scalar LTV system: ẋ1 = p ′ 3 2 b c (t) a (t) c (t) x + u; 1 2/3 b′ a (t) y = x1 (32) ′ Here a, b , c are the physical parameters which may represent quantities like mass, friction coefficient, resistance, ′ capacitance etc. a and c are time varying while b is an uncer ∗ tain parameter. a and c varies as: a (t) = a∗ − a2 exp − τta ,   ∗ ∗ ′ c (t) = c2 + c2 exp − τtc and b lies in the uncertainty h i ∗ set b∗ , 3b2 . Here physical parameter Θ = [a, b, c] and hence p = 3. Fromh the time i variation hof∗ a and i c we can ∗ conclude that a ∈ a2 , a∗ and c ∈ c2 , c∗ . Therefore, 14 Semi-Infinite 15 For optimization problems are the ones with infinite constraints. a given L, P is not unique. (a∗ , 3b∗ /2, c∗ ) ! ′ a, b , c ! " ′ (a1 , b) = F a, b , c " (a∗ /2, b∗ , c∗ /2) (a1 , b) b) a) c) Figure 8. a) Uncertainty set S associated with the physical parameters. The coordinates of two diagonally opposite vertices is shown in the figure. b) Uncertainty set L associated with the LTV system characterized by equation (15). c) Bounding Convex Polytope P of set L. The red lines are the edges of the polytope P. P has 5 vertices which forms the elements of set V. h ∗ i h ∗ i h i ∗ S ∈ a2 , a∗ × b∗ , 3b2 × c2 , c∗ which is an hypercube as shown in Fig. 8a. To find set L first note that for this example N = 1 as the system described by equation (32) is scalar. Therefore T T L ⊂ R2 and  every element [a1 , b] ∈ L is a map [a1 , b] = ′ F a, b , c where " a3 c2 #   ′ ′ F a, b , c = √b ′ (33) b c a2/3 One method to obtain set L is to divide the hypercube S into uniform grids and then map each grid using equation (33). Set L for this example is shown in Fig. 8b. Though this method to obtain L from S is computationally expensive, one must appreciate that for a general system there is no other elegant method. Now we want to find a convex polytope P which bounds L. One such polytope is shown in Fig. 8c. This polytope has m = 5 vertices. In practise polytope P can be found using convex-hull functions available widely in many programming languages. One popular reference is [29]. We now discuss how to calculate δA and δB for this example. At this point one must appreciate that in most practical scenario the controller designer may not explicitly know the equations governing the rate of change of physical parameters, i.e. say in this example the designer may not know a (t) and c (t) explicitly. However it seems practical to assume knowledge of the bounds on the rate of change of physical parameters, i.e. controller designer may know that a∗ c∗ 0 < ȧ ≤ 2τ and − 2τ ≤ ċ < 0. There is no standard way to a c calculate δA and δB from these bounds. Also dependent on the method used, one may get different estimate of δA and δB. For this example δA and δB is calculated as follows.   2 2   3 2  |3c a ȧ+2a3 cċ| a(t) c(t) d δA = max dt = max b′ b′ = max(3 sup(c)2 sup(a)2 max(ȧ), 2 sup(a)3 sup(c) max(|ċ|)) inf (b′ ) ∗ = = ∗ a c max(3(c∗ )2 (a∗ )2 ( 2τ ), 2(a∗ )3 (c∗ )( 2τ )) a c (b∗ )2 (a∗ )3 (c∗ )2 2(b∗ )2 max  3 2 τa , τc  10 δB = max  √ d dt b′ c(t) a(t) √ 2/3 √   ċ b′ − 23a5c/3ȧ + 2√ca 2/3  √ ∗ a∗ p 2 sup(c)( 2τ (c ) ) a = sup (b′ ) + √ 2τc 5/3 = max 3 inf(a) √ c∗  m X 2 inf(c) inf(a) 2/3 i=1 m X  2 1 2/3 3τc + 25/6 τc (a∗ ) Lemma 3: Under Assumption 4, if for a given P ∈ RN ×N and s ∈ R there exist a Ki ∈ [km , kM ] s.t. = T (Ai − Bi Ki C) P + P (Ai − Bi Ki C)  sI, (Ai , Bi ) ∈ V then there exist a KAB ∈ [km , kM ] s.t. T (A − BKAB C) P +P (A − BKAB C)  sI ∀ (A, B) ∈ L Proof: We first define a N × N matrix Γ all elements of which are 0 except its (N, 1) element which is 1. Also note that all (A, B) ∈ L can be written as a convex combination of the elements V of the convex polytope P. Mathematically, ∀ (A, B) ∈ L there exists scalars θi ≥ 0, i = 1, 2, . . . m s.t. (A, B) = m X θi (Ai , Bi ) and m X θi = 1 i=1 i=1 where, (Ai , Bi ) ∈ V ∀i = 1, 2, . . . , m. Now, T (A − BKAB C) P + P (A − BKAB C)

= AT P + P A − bKAB ΓT P + P Γ = m X θi ATi P i=1

T + P Ai − Γ P + P Γ KAB m X θi bi (34) i=1 Equation (34) is possible because BC = b (t) Γ owing to the special structure of B and C matrix. Using the inequality T (Ai − Bi Ki C) P + P (Ai − Bi Ki C)  sI, (Ai , Bi ) ∈ V in equation (34) we have T (A − BKAB C) P + P (A − BKAB C) # "m m X X
θi bi ΓT P + P Γ + sI θi Ki bi − KAB  (35) i=1 i=1 Choosing, KAB = m X θi K i b i i=1 m X in inequality (35) yields θi b i i=1 ≤ KAB ≤ m X θ i k M bi i=1 m X θ i bi ⇒ km ≤ KAB ≤ kM i=1 This concludes the proof of Lemma 3. The importance of Lemma 3 is that it reduces the semi-infinite optimization problem posed in Problem 1 into a finite optimization problem. This results into the most important result of this section. Theorem 2: [C1] holds if the optimal value of the problem minimize: s subject to: T (Ai − Bi Ki C) P + P (Ai − Bi Ki C)  sI ∀ (Ai , Bi ) ∈ V km ≤ Ki ≤ kM , i = 1, 2, . . . , m P ≻0 with design variables s ∈ R, P ∈ RN ×N and Ki ∈ R, i = 1, 2, . . . , m is negative. Proof: The proof follows from Problem 1 and Lemma 3. Note that while Problem 1 is a necessary and sufficient condition, Theorem 2 is a sufficient condition. This is due to the fact that L ⊆ P leading to some conservativeness in the optimization problem proposed in Theorem 2. Theorem 2 poses the classical BMI Eigenvalue Minimization Problem (BMI-EMP) in variables P and Ki . As such BMI’s are non-convex in nature leading to multiple local minimas. Several algorithms to find the local minima exist in literature (see [30], [31]). Algorithms for global minimization of BMI’s is rather rare and have received attention in works like [32], [33]. Our approach is similar to [32] which is basically a Branch and Bound algorithm. Such an algorithm works by bounding s by a lower bound ΦL and an upper bound ΦU , i.e. ΦL ≤ s ≤ ΦU . The algorithm then progressively refines the search to reduce ΦU − ΦL . Our main algorithm consists of Algorithm 4.1 of [32] (Page 4). The Alternating SDP method mentioned in [34] (Page 2) is used for calculating ΦU . For calculating ΦL we have used the convex relaxation technique first introduced in [33] (also discussed in [34], equation (9)). In Appendix B we present a working knowledge of our algorithm. For detailed explanation the readers may refer the corresponding work [32], [34]. Theorem 2 poses an optimization problem with (Ai , Bi ) ∈ V, km and kM as inputs and P , s and Ki , i = 1, 2, . . . , m as outputs. But km and kM are not known. An initial esRH RH is obtained by and kM = kM timate of km = km using Routh-Hurwitz criteria17 for each (Ai , Bi ) ∈ V. Let Ki = KiRH , i = 1, 2, . . . , m be the output of the optimization
∗ = min KiRH problem with this initial estimate. Let km 1≤i≤m
∗ = max KiRH , then the following holds: and kM 1≤i≤m T (A − BKAB C) P + P (A − BKAB C)  sI Now all we need to do is to prove that the chosen KAB lies in the interval [km , kM ]. We know that km ≤ Ki ≤ kM . Therefore under Assumption 416 , 16 Without θ i bi i=1  5/3 θ i k m bi Assumption 4, the denominator Pm i=1 θi bi may become 0. RH RH ∗ RH ∗ ≤ . This is because km ≤ kM and kM ≥ km 1) km RH KiRH ≤ kM , ∀ i = 1, 2, . . . , m and hence

RH RH ≤ min KiRH ≤ max KiRH ≤ kM km 1≤i≤m 1≤i≤m 17 Routh-Hurwitz criteria is used to find the bounds on the feedback gains for which a SISO LTI system is closed loop stable. Refer [35] for details. 11 2) The outputs, P , s and Ki , obtained by running the opRH RH timization algorithm with a) km = km and kM = kM ∗ ∗ or b) km = km and kM = kM , will be the same. This is because the gains KiRH obtained by running the optimization RH RH algorithm with km = km and kM = kM also satisfies the ∗ RH ∗ bounds km ≤ Ki ≤ kM . ∗ ∗ Therefore if we choose km = km and kM = kM ∗ ∗ we would get
a smaller RGS gain set, i.e. (kM − km ) ≤ RH RH kM − km , without compromising the convergence coefficient s. A smaller RGS gain set will ease the controller design in an analog setting. We now give a bound on λL M (defined in Theorem 1). It is not possible to calculate λL M with the formula given in L Theorem 1 as it will involve search over the dense set SQ . Define a setn T P SQ = Q : Q = (A − BKC) P + P (A − BKC) o (A, B) ∈ P, K ∈ [km , kM ] L Let λP M = max λM (Q). As L ⊆ P it is obvious that, SQ ⊆ P Q∈SQ P L P SQ (SQ defined in Theorem 1). Therefore18 , λL M ≤ λM . Thus L P λM gives an estimate of λM by upper bounding it. It can be shown that for a scalar gain K and the specific structure of P B and C (refer equation (15)), it can be proved that SQ is compact convex set. Also λM (Q) is a convex function for all Q ∈ S+ (refer Page 82 of [36]). It is well known that global maxima of a convex function over a convex compact set only occurs at some extreme points of the set (refer [37]). Thus P the problem of maximizing λM (·) over Q ∈ SQ reduces to V maximizingnλM (·) over Q ∈ SQ where T V SQ = Q : Q = (A − BKC) P + P (A − BKC) o (A, B) ∈ V, K ∈ {km , kM } P . This leads to the following formula the set of vertices of SQ P λL M ≤ λM = max λM (Q) V Q∈SQ (36) Inequality (36) can be used to obtain an estimate of λL M. Remark 8: As mentioned in the beginning of Section IV, for a controller to be implementable in analog framework it has to be simple. Though the synthesis of RGS parameters (km , kM , T , α, γ and P ) is complex, one must appreciate that designing a controller is a ’one time deal’. RGS in itself is a simple gain-scheduled controller governed by equation (16). V. E XAMPLE Parallel Plate Electrostatic Actuator (PPA) shown in Fig. 9a forms a vital component of several miniaturized systems. We perform regulatory control of PPA to show effectiveness of the proposed analog architecture and RGS strategy. PPA’s as described in [20] (page 183) follows the following dynamics mÿ + bẏ + ky = εA 2 (G − y) 2 2 Vs (37) 18 This is more like saying that the maximum of a function (λ M (·) here) P here) will be greater than the maximum over a smaller over a bigger set (SQ L here). set (SQ which is nonlinear in nature. Plant parameter includes spring constant κ, damping coefficient b, moving plate mass and area m and A respectively, permittivity ε and maximum plate gap G. As we are interested in only regulatory control, we use the following linearized model # #  "   " 0 0 1 x1 x˙1 √ + = 2εAκGo ur (38) b o) x˙2 x2 − κ(G−3G −m m(G−Go ) m(G−Go ) where, x1 is the displacement from the operating point Go . Plant output is the moving plate position y = Go + x1 . Plant input is Vs = Vb + Vu . Comparing with RGS theory, Vb = ub , the bias voltage to maintain Go as the operating point and Vu = ur = KVe = K (Go − y) = −Kx1 , the regulation voltage supplied by the Memristive AGC. Plant parameter includes spring constant κ, damping coefficient b, moving plate mass and area m and A respectively, permittivity ε and maximum plate gap G. For Go > G 3 , the system has an unstable pole. We perform regulation around Go = 2G 3 . The true plant parameters are, m = 3×10−3 kg, b = 1.79× −3 10−2 N sm−1 , G = 10−3 m, A = 1.6×10 m2 . G and A are  uncertain but G ∈ SG = 0.5 2.0 × 10−3 m,  lie in the  set−3 A ∈ SA = 1.2 1.8 × 10 m2 . ε varies due to surrounding condition as ε (t) = 5εo + 1.5εo sin (7.854t) where εo is the permittivity of vacuum. Spring loosening causes κ to decrease as κ (t) = 0.08+0.087e−0.8t N m−1 . Now we will discuss the steps involved in implementing RGS in an analog framework. Step 1 (Identify S): S is the uncertainty set of physical param-  3.5ε 6.5εo eters first defined in Page 8. Define set S = o ε   and Sκ = 0.08 0.167 N m−1 . Then the ordered pair (G, A, ε, κ) ∈ S = SG × SA × Sε × Sκ . Note that here p = 4. Step 2 (Find P): To do this we numerically map S to L (as shown in Example 2) and then use convhulln function of MATLAB to find a convex polytope P s.t. L ⊆ P. In this case P consist of m = 4 vertices. We explicitly don’t mention the computed P for the sake of neatness. Step 3 (Compute P , s, α, km , kM ): Solving optimization problem proposed in Theorem 2 we get, s = 0.917, km = 0.9937 0.0757 8600 and kM = 86000 and P = . Here, 0.0757 0.0895 λM (P ) = 1, hence α = λMs(P ) = 0.917. Step 4 (Compute λm (P ), λL M , δ): For the calculated P , λm (P ) = 0.083. From equation (36), λL M = 29.1. We will now calculate δA and δB for this example. As mentioned in Example 2 the controller designer does not know κ (t) and ε (t) explicitly but they do know the bounds on κ̇ and ε̇ which for this example is: (−0.087 × 0.8) ≤ κ̇ < 0 and − (1.5 × 7.854 × εo ) ≤ ε̇ ≤ (1.5 × 7.854 × εo ). First note that (2 × 2) element of the system matrix of the linearized PPA model (described by equation (38)) is not time varying. Hence δA can be simply written as    d o) − κ(t)(G−3G δA = max dt m(G−Go ) = max =  d dt 3 max(|κ̇|) m  − κ(t)(G−3( 2G 3 )) m(G− 2G 3 )  12 100 b 50 G m A + y x1 (in µm) κ Operating Point Go x1 Moving Plate Vs 0 −50 −100 Fixed Plate −150 a) −200 x1 + K km P22 b) 0 R1 + Vsh R2 8 TC1 − Ė E Ė V + TC3 Go R − VT + y Ve VC Memristive Analog Gain Controller Vu VL1 VDD −VDD 3 3.5 0 −8 0.5 x 10 1 1.5 2 2.5 e) t (in s) 3 3.5 − E + Voltage Toggler + 2.5 6 Hysteresis Block R 2 d) t (in s) Scan Mode Rest Mode 2 R3 1.5 4 T1 −VDD TC2 1 6 x2 d dt 0.5 8 E 2P12 E 0 10 P11 4 2 VL2 0 0 0.5 1 c) 1.5 f) 2 2.5 t (in s) 3 3.5 Figure 9. a) Schematic of PPA. b) Analog Implementation of E = xT P x for P ∈ R2×2 . Note that, P12 = P21 as P ∈ S. c) Analog Implementation of RGS. d), e), f) Plots of plate position error x1 (from operating point Go ), normalized RGS gain kK and Lyapunov Function E = xT P x with respect m to time. δB = max = max = 1 m = 1 m = 1 m q q q  d dt d dt √ q 2ε(t)Aκ(t)Go m(G−Go ) 2ε(t)Aκ(t)( 2G 3 ) m(G− 2G 3 ) 12 sup(SA ) inf(SG ) max 12 sup(SA ) inf(SG ) max  12 sup(SA ) inf(SG )  !!  √ ˙ εK   κ√ ε̇+ε √κ̇ 2 ε κ sup(Sκ ) max(|ε̇|)+sup(Sε ) max(|κ̇|) 2 √ inf(Sε ) √ inf(Sκ )  Substituting the values in the above equation of δA and δB we get δ = δA + kM δB = 1351. Step 5 (Compute T and γ): We arbitarily choose γ = 0.5. Substituting γ = 0.5 in inequality (29) we get T < 1.66 × 10−7 s. Hence we choose T = 10−7 s. Step 6 (Analog Design): Fig. 9b and 9c combined shows the analog implementation of RGS. The error voltage Ve is the plate position error, i.e. Ve = −x1 = (Go − y). The gain control voltage VC is controlled by the Hysteresis Block and the Voltage Toggler block. In Rest Mode VC = 0 thereby ensuring that the gain M is constant (refer equation (12)). In Scan Mode VC ∈ {−VDD , VDD }, VC = VDD to scan from S S S S Rof f to Ron and VC = −VDD to scan from Ron to Rof f (refer equation (12)). RC controls Ṁ, the rate of change of DD T gain (refer equation (12)). Setting RC = VQ ensures a S M scan time of T . The derivation of RC is simple. Note that the S S resistance of memristor will change from Rof f to Ron if we pass a charge QSM through it. We want this change to happen QS in time T and hence the desired current is TM . As VC = VDD DD T . = VQ in scan mode, the required resistance RC = QVSDD S ( M/T ) M The Hysteresis Block is a conventional inverting schmitt 2(VDD −α) DD −α) trigger. Tuning R1 = 2(V α(1−γ) R2 and R3 = α(1+γ) R2 19 ensures that the schmitt trigger’s output goes from VDD to Ė −VDD at Ė E = −γα and from −VDD to VDD at E = −α. Due to the transistor arrangement in Hysteresis Block: 1) VC = 0 when Vsh = VDD . Therefore Vsh = VDD implies Rest Mode. 2) VC = VT when Vsh = −VDD . It will be explained next that VT ∈ {−VDD , VDD }. Therefore Vsh = −VDD implies Scan Mode. So we can conclude that the Hysteresis Block Ė = −γα and goes from Rest Mode to Scan Mode when E Ė from Scan Mode to Rest Mode when E = −α. This is in accordance with the hysteresis shown in Fig. 6a. VT , the output of the Voltage Toggler block, toggles between −VDD and VDD . Recalling equation (12), this will result in the  memristive gain control block reflect between  gain of the S S Rof Ron f . We now explain the working of this block. Say RI , RI VT = VDD and the zone indicating voltages (refer Fig. 3) VL1 = VL2 = VDD . As VL1 = VL2 = VDD , transistors TC2 and TC3 are OF F . The voltage V + = V2T = VDD 2 . Since V + = VDD > 0, V = V . This shows that the output of T DD 2 Voltage Toggler block is stable. As VT = VDD , memristor’s S resistance M will decrease till M = Ron (see equation (12)) at which point VL1 = −VDD and VL2 = VDD (refer Case 2 of III-D). TC3 will be momentarily ON making V + = −VDD and hence driving VT to −VDD . Then M will increase from S Ron , making VL1 = VL2 = VDD and hence driving TC3 to + = OF F state. When this occurs V + = V2T = − VDD 2 . As V VDD − 2 < 0, VT = −VDD . Similar momentary transition of TC2 to ON state will toggle VT from −VDD to VDD when S M = Rof f. Several plots corresponding the regulatory performance of PPA under RGS control strategy is shown in Fig. 9d, e, f. It is interesting to observe that an LTV system can have multiple rest-scan cycle (see Fig. 9e). This is because for a time varying system a stabilizing gain at a given instant may not be the stabilizing gain at a later instant due to the change in system dynamics. Unlike LTV system, a LTI system will have only 1 rest-scan cycle. Remark 9: In this example T is very low which may seem to challenge analog design. However for all practical purposes it is not so. For the sake of simulation we choose a fictitious time variation of κ (t) and ε (t) which is quite fast compared to that found in nature. Therefore δA and δB is high (refer Step 4) resulting in a high δ and hence a low scan time T (refer inequality (29). In practice, time variation of a system caused by ageing effect and atmospheric variation is a slow process. Hence, T will be much higher. Remark 10: To control an array of miniaturized devices (say PPA) one can reduce the circuitry required by identifying the components of the circuit which can be common for the entire array. For example, Synchronization Block can be common for the entire array. Synchronization of each pair of coupled Memristor Gain Block and Integrator (refer Fig. 3) can be done in a time multiplexed manner, i.e. each pair of coupled Memristor Gain Block and Integrator is synchronized using 19 The designer is free to choose the resistance R2 . 13 one Synchronization Block. The oscillator shown in the circuit of Fig. 3 can also be common for the entire array. In equation (40) δA = Ȧ (t∗ ) δt and δB = Ḃ (t∗ ) δt are infinitesimal change in A and B respectively. Note that ∆K is a scalar. Hence, substituting equation (39) in (40) yields T VI. D ISCUSSION To the best of authors knowledge this paper is one of the first attempts towards understanding the use of memristor in control applications. A memristive variable gain architecture has been proposed. We then propose (and prove) RGS control strategy which can be implemented using this framework. Simplicity of RGS control strategy is demonstrated using an example. The extension of this work can take two course. From Circuit Standpoint one may try to design an analog circuit which mimics the circuit shown in Fig. 2 but with a lesser number of op-amps. Since the synthesis of memristor is still an engineering challenge, one may speculate regarding the use of variable CMOS resistors (refer [21]) to implement the analog gain controller proposed in Section III. Two milestones have to be addressed before RGS is practically feasible: 1) RGS needs information about the states x which is obtained by differentiating the output y N − 1 times. But differentiation might amplify the noise. 2) RGS relies on the knowledge of E which is obtained by performing xT P x using analog circuitry. Such analog implementation will be easier if P is sparse. Hence from Control Theoretic Standpoint, addressing these two issues will be the first target of the authors. Later point has been addressed in [38]. Extending RGS to SISO non-linear and in general MIMO systems would be the next step. It would also be interesting to explore other simple control strategies (like [39]) which can be implemented in analog framework. A PPENDIX A P ROOF OF L EMMA 1: D RIFTING NATURE OF S TABILIZING G AIN S ET KP,s (t) To prove Lemma 1 we will take the following steps: 1) Pick a gain K ∈ KP,s (t). 2) Prove that if δt → 0 then there exist a ∆K → 0 s.t. (K + ∆K) ∈ KP,s (t + δt). 3) As ∆K → 0, (K + ∆K) ∈ KP,s (t). Hence the gain K + ∆K belongs to both the sets, KP,s (t) T and KP,s (t + δt). This implies that KP,s (t) KP,s (t + δt) 6= ∅. Now we proceed with the proof. We first pick a gain K ∈ KP,s (t). As KP,s (t) 6= ∅, ∀t ≥ 0 such a gain will exist. For a time t = t∗ there exists a gain K ∗ ∈ KP,s (t∗ ) and a scalar s∗ > s s.t. the following equality holds T (A∗ − B ∗ K ∗ C) P + P (A∗ − B ∗ K ∗ C) = −s∗ I (39) In equation (39) A∗ = A (t∗ ) and B ∗ = B (t∗ ). Equation (38) directly follows from the very definition of stabilizing gain set (defined in page 7). Lets say that at time t = t∗ + δt there exist a K ∗ + ∆K s.t. T [(A∗ + δA) − (B ∗ + δB) (K ∗ + ∆K) C] P +P [(A∗ + δA) − (B ∗ + δB) (K ∗ + ∆K) C] = −s∗ I (40) ∆KR = (δA − δBK ∗ C) P + P (δA − δBK ∗ C) (41) T where, R = C T (B ∗ + δB) P + P (B ∗ + δB) C. Hence if ∆K satisfies equation (41) then K ∗ + ∆K ∈ KP,s (t∗ + δt). Now all we need to do is to prove that ∆K is infinitesimal, i.e. ∆K → 0 if δA → 0 and δB → 0. Taking norm on both side of equation (41) we get, ∆K kRk ≤ ∆K ≤ 2 kP k (kδAk + K ∗ kδBk kCk) 2 kP k (δA + K ∗ δB ) δt kRk (42) Since kRk is finite it is obvious from inequality (42) that ∆K → 0 as δt → 0. This concludes the proof. A PPENDIX B G LOBAL S OLUTION OF BMI-EMP Theorem 2 involves solving the following BMI-EMP optimization problem OP1: minimize: s subject to: T (Ai − Bi Ki C) P + P (Ai − Bi Ki C)  sI ∀ (Ai , Bi ) ∈ V km ≤ Ki ≤ kM , i = 1, 2, . . . , m P ≻0 We will first justify why OP1 is called BMI-EMP. Consider a matrix inequality of type T (A − BKC) P + P (A − BKC)  sI (43) Given a set of A, B, C, K and P , the least possible value of s which will satisfy inequality (43) is indeed the largest T eigen value of the matrix (A − BKC) P + P (A − BKC). So if we exclude the inequalities km ≤ Ki ≤ kM and P ≻ 0, then OP1 can be equaly casted as   T min max λM (Ai − Bi Ki C) P + P (Ai − Bi Ki C) P, Ki (Ai , Bi )∈V which is basically a Largest Eigenvalue Minimization Problem or just “EMP”. Now consider the function   T Λ (P , K) = λM (A − BKC) P + P (A − BKC) T The matrix Q (P, K) := (A − BKC) P +P (A − BKC) is Bilinear in the sense that it is linear in K if P is fixed and linear in P if K is fixed. As Q ∈ S, the function Λ (P , K) = λM (Q (P, K)) is convex in K if P is fixed and convex in P if K is fixed. We are interested in the global solution of BMI-EMP as a smaller value s will ensure better convergence of the LTV/LTI system. In the following we will provide a sketch of the work done in [32], [34], [33] which will be just enough to design an algorithm for global solution of BMI-EMP. However we don’t provide detailed explanation of the algorithm for which the reader may refer [32]. Before proceeding forward we would like to cast OP1 in a form which can be handled by a numerical solver. Observe that 14 the third constrain of OP1 is a strict inequality which demands that the Lyapunov Matrix P has to be positive definite (not positive semi-definite). Such a strict inequality will impose numerical challenge and hence we replace it with the nonstrict inequality P  µp I where, 0 < µp ≪ 1. Without any loss of generality: 1) We constrain the kP k ≤ 1 by imposing the constrain P  I. 2) We normalize the RGS gain set. We get the following optimization problem, OP2: minimize: s subject to: T ATi P +P Ai −(Bi C) Ki P −Ki P (Bi C)  sI ∀ (Ai , Bi ) ∈ V µk ≤ Ki ≤ 1, i = 1, 2, . . . , m µp I  P  I km In the above  problem, µk =  kM < 1. We also redefine C as C := kM · · · 0 0 , to neutralize the effect of normalizing RGS gain set.  T We now define two vectors: P = p1 p2 · · · pnp containing the np = N (N2+1) distinct elements of symmetric  T matrix P and K = K1 K2 · · · Km the m normalized RGS Gains of OP2. We also define two sets XP and XK as follows n m XP := [−1, 1] p ; XK := [µk , 1] Note that XP is a np dimensional unit hypercube such that20 P ∈ XP and XK is a m dimensional hyper-rectangle such that K ∈ XK . We also define smaller hyper-rectangle’s QP ⊆ XP , QK ⊆ XK and Q ⊆ XP × XK as follows   n    n  QP := L1P , UP1 × L2P , UP2 × . . . × LPp , UP p    2  m 2 1 × . . . × [Lm × LK , U K QK := L1K , UK K , UK ] Q := QP × QK i Obviously −1 ≤ LiP ≤ UPi ≤ 1 and µk ≤ LiK ≤ UK ≤ 1. Consider the following “constrained version” of OP2 where (P , K) is only defined in the small hyper-rectangle Q. OP3: minimize: s subject to: T ATi P +P Ai −(Bi C) Ki P −Ki P (Bi C)  sI ∀ (Ai , Bi ) ∈ V µk ≤ Ki ≤ 1, i = 1, 2, . . . , m µp I  P  I (P , K) ∈ Q The input to OP3 is the set Q and its output is s∗ (Q) which is a function of Q. We want to bound s∗ (Q) by an upper and a lower bound as follows: ΦL (Q) ≤ s∗ (Q) ≤ ΦU (Q) (44) The convex relaxation technique introduced in [33] (also discussed in [34], equation (9)) has been used to get the lower bound ΦL (Q). We replace the nonlinearity Ki P with a new matrix Wi . As Wi is a symmetric matrix matrix of order N we can represent it by a vecT w1i w2i · · · wnp i containing the np tor Wi = 20 All the element of the Lyapunov Matrix P will be in the range [−1, 1] as P  I according to OP2. Algorithm 1 Alternating SDP Method

1. Set K(0) =Centroid of QK . s(0) , P(0) = OP3 Q, K(0) .
2. if s(0) = ∞ 3. return ∞. 4. else 5. Set δ > 0 and k = 0. 6. do {

7. s(k+1) , P(k+1)
= OP3 Q , K(k)
8. s(k+1) , K(k+1) = OP3 Q , P(k+1) 9. k = k + 1.
10. } while s(k−1) − s(k) < δ s(k) 11. return s(k) . 12. end elements of Wi . We T now define a matrix W = distinct T W1T W2T · · · Wm . If OP3 is expanded in terms of pj and Ki then the element wji of W is constrained by the equality wji = Ki pj . Rather than imposing this equality constrain we let wji to be a free variable which can take any value in the set W (Q) defined as   (P , K) ∈ Q        wji ≥ LiK pj + LjP Ki − LiK LjP    i i pj + UPj Ki − UK UPj W (Q) := W wji ≥ UK   i i    LjP  wji ≤ UK pj + LjP Ki − UK     j j i i wji ≤ LK pj + UP Ki − LK UP By performing the convex relaxation stated above we get the following optimization problem. OP4: minimize: s subject to: T ATi P + P Ai − (Bi C) Wi − Wi (Bi C)  sI ∀ (Ai , Bi ) ∈ V µp I  P  I µp Ki Iµk  Wi  Ki I, i = 1, 2, . . . , m (P , K) ∈ Q W ∈ W (Q) The input to OP4 is the set Q and its output is ΦL (Q). As OP4 is a relaxed version of OP3, ΦL (Q) ≤ s∗ (Q). Note that OP4 is a convex optimization problem, more specifically a Semi-Definite Program(SDP) which can be solved by numerical solvers like CVX[40]. We now concentrate on defining the upper bound ΦU (Q). Any local minima of OP3 can indeed be the upper bound ΦU (Q). We use the Alternating SDP Method discussed in [33], [41] . Alternating SDP method relies on the Bi-Convex nature of OP3, i.e. OP3 becomes a convex problem (more specifically SDP) in P with K fixed or a convex problem in K with P fixed. Alternating SDP Method is summarized in Algorithm 1. We represent by OP3 (Q , K′ ) the opti′ mization problem obtained by fixing K = K in OP3 and ′ OP3 (Q , P ) the optimization problem obtained by fixing ′ P = P in OP3. The input to OP3 (Q , K′ ) is the small ′ hyper-rectangle Q and the fixed RGS gain K while the input ′ to OP3 (Q , P ) is the small hyper-rectangle Q and the fixed ′ Lyapunov Matrix P . The outputs of OP3 (Q , K′ ) are the optimized s and P while the outputs of OP3 (Q , P′ ) are the optimized s and K. 15 Algorithm 2 Branch and Bound 1. Set ǫ > 0 and k = 0. 2. Set Q0 = XP × XK and G0 = {Q0 }. 3. Set L0 = ΦL (Q0 ) and U0 = ΦU (Q0 ). 4. while (Uk − Lk < ǫ)
5. Select Q̄ from Gk such that Lk = ΦL Q̄ . 6. Set Gk+1 = Gk − Q̄ . 7. Split Q̄ along its longest egde into Q̄1 and Q̄2 . 8. f or (i = 1, 2)

9. if ΦL Q̄i ≤ Uk
10. Compute ΦU Q̄ Si . Q̄i . 11. Set Gk+1 = Gk 12. end 13. end 14. Set Uk+1 = min ΦU (Q). 15. 16. Q∈Gk+1 P runing : Gk+1 = Gk+1 − {Q : ΦL (Q) > Uk+1 }. Set Lk+1 = min ΦL (Q). Q∈Gk+1 17. Set k = k + 1. 18. end Under the various definations introduced above the Branch and Bound Algorithm [32] calculate the global minima of OP2 to an absolute accuracy of ǫ > 0 within finite time. The psuedocode of Branch and Bound method is given in Algorithm 2. A PPENDIX C C ONTROL T HEORETIC D EFINITIONS AND C ONCEPTS The control theoretic analysis presented in this paper relies on definitions and theorems from three broad areas. This concepts are standard and can be easily found in books like [42], [36]. A. Real Analysis Defintion 1 (Cartesian Product of Sets): Cartesian Products of two sets A and B, denoted by A × B, is defined as the set of all ordered pairs (a, b) where a ∈ A and b ∈ B. Notion 1 (Connected Set): A set is said to be connected if for any two points on that set, there exist at-least one path joining those two points which also lies on the set. Note that this is not the "general definition" of connected set. Notion 2 (Compact Set): In Eucledian Space Rn , a closed and bounded set is called a compact set. Also a compact set in Rn is always closed and bounded. Note that this is not the "general definition" of compact set. B. Linear Algebra Definition 2 (Eucledian Norm of a Vector): For a √vector x ∈ Rn , eucledian norm kxk is defined as kxk := xT x. Throughout the entire paper “norm” means “eucledian norm” unless mentioned otherwise. Definition 3 (Induced Spectral Norm of a Matrix): For a matrix A ∈ Rm×n , Induced Spectral Norm kAk is defined as kAk := kAxk x∈Rn −{0} kxk sup It can equally be defined as kAk := sup kAxk kxk=1 Definition 4 (Positive (Negative) Definite (Semi-Definite) Matrix): A square matrix A ∈ Rn×n is said to be positive definite if xT Ax > 0, ∀ x ∈ Rn −{0} and is said to be positive semi-definite if xT Ax ≥ 0, ∀ x ∈ Rn . A matrix A is said to be negative definite is −A is positive definite and is said to be negative semi-definite is −A is positive semi-definite. Linear Algebra Theorems: 1) Properties of norms: a) For a vector x ∈ Rn , if kxk = 0 then x = 0. Also if kxk → 0 then x → 0. b) For any two matrix A, B ∈ Rn×m , kA + Bk ≤ kAk + kBk. c) For any matrix A ∈ Rn×m and a scalar α, kαAk = α kAk. d) For any two matrix A, B ∈ Rn×m , kABk ≤ kAk kBk. e) For any matrix A ∈ Rn×m and a vector x ∈ Rm , kAxk ≤ kAk kxk. 2) Eigenvalues of symmetric positive (negative) definite matrix are positive (negative). 3) For any symmetric matrix A ∈ Rn×n and a vector x ∈ Rn , 2 λm (A) kxk ≤ xT Ax ≤ λM (A) kxk 2 4) For any symmetric positive definite (semi-definite) matrix A, kAk = λM (A). 2 5) For any square matrix A, xT Ax ≤ kAk kxk C. Convex Analysis Definition 5 (Convex Set): A set A is said to be convex if for any x , y ∈ A, θx + (1 − θ) y ∈ A , ∀ θ ∈ [0 , 1]. Definition 6 (Convex Function): Let A ∈ Rn be a convex set. A function f : A → R is convex if f (θx + (1 − θ) y) ≤ θf (x) + (1 − θ) f (y) for all x , y ∈ A and for all 0 ≤ θ ≤ 1. Definition  7 (Convex Combination and Convex Hull): Given a set A = a1 a2 · · · an , its convex combination are those elements which can be expressed as θ1 a1 + θ2 a 2 + · · · + θn an where θ1 + θ2 + · · · + θn = 1 and θi ≥ 0, i = 1, 2, . . . , n. The set of all convex combination of set A is called the convex hull of A. Convex Hull of set A is indeed the smallest convex set containing A. Definition 8 (Semi-Definite Program): A semi-definite program or SDP is an optimization problem in variable x =  T x1 x2 · · · xn ∈ Rn with the generic structure minimize : cT x subject to : F0 + x1 F1 + x2 F2 + · · · + xn Fn  0 Ax = b where F0 , F1 , F2 , . . . , Fn ∈ S, c ∈ Rn , A ∈ Rp×n and b ∈ Rp . 16 ACKNOWLEDGMENT We would like to thank Prof. Radha Krishna Ganti for providing invaluable help in areas related to convex optimization. R EFERENCES [1] L. Chua, “Memristor-the missing circuit element,” IEEE Trans. Circuit Theory, vol. 18, no. 5, pp. 507–519, 1971. [2] D. B. Strukov, G. S. Snider, D. R. Stewart, and R. S. Williams, “The missing memristor found,” Nature, vol. 453, no. 7191, pp. 80–83, 2008. [3] Y. Ho, G. M. Huang, and P. Li, “Dynamical properties and design analysis for nonvolatile memristor memories,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 58, no. 4, pp. 724–736, 2011. [4] Y. N. Joglekar and S. J. Wolf, “The elusive memristor: properties of basic electrical circuits,” Eur. J. Phys., vol. 30, no. 4, pp. 661–675, 2009. [5] F. Corinto and A. Ascoli, “A boundary condition-based approach to the modeling of memristor nanostructures,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 60, pp. 2713–2726, 2012. [6] S. Kvatinsky, E. G. Friedman, A. Kolodny, and U. C. Weiser, “Team: Threshold adaptive memristor model,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 60, no. 1, pp. 211–221, 2013. [7] S. H. Jo, T. Chang, I. Ebong, B. B. Bhadviya, P. Mazumder, and W. Lu, “Nanoscale memristor device as synapse in neuromorphic systems,” Nano letters, vol. 10, no. 4, pp. 1297–1301, 2010. [8] Y. V. Pershin and M. Di Ventra, “Experimental demonstration of associative memory with memristive neural networks,” Neural Networks, vol. 23, no. 7, pp. 881–886, 2010. [9] A. Sheri, H. Hwang, M. Jeon, and B. Lee, “Neuromorphic character recognition system with two pcmo-memristors as a synapse,” IEEE Trans. Ind. Electron., vol. 61, pp. 2933–2941, 2014. [10] Y. V. Pershin, S. La Fontaine, and M. Di Ventra, “Memristive model of amoeba learning,” Phys. Rev. E, vol. 80, no. 2, p. 021926, 2009. [11] G. Z. Cohen, Y. V. Pershin, and M. Di Ventra, “Second and higher harmonics generation with memristive systems,” Appl. Phys. Lett., vol. 100, no. 13, p. 133109, 2012. [12] Y. V. Pershin and M. Di Ventra, “Practical approach to programmable analog circuits with memristors,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 57, no. 8, pp. 1857–1864, 2010. [13] S. Shin, K. Kim, and S. Kang, “Memristor applications for programmable analog ics,” IEEE Trans. Nanotechnol., vol. 10, no. 2, pp. 266–274, 2011. [14] D. Jeltsema and A. Doria-Cerezo, “Port-hamiltonian formulation of systems with memory,” Proc. IEEE, vol. 100, pp. 1928–1937, 2012. [15] D. Jeltsema and A. J. van der Schaft, “Memristive port-hamiltonian systems,” Math. Comput. Model. Dyn. Syst., vol. 16, pp. 75–93, 2010. [16] J. Hu and J. Wang, “Global uniform asymptotic stability of memristorbased recurrent neural networks with time delays,” in Proc. Int. Joint Conf. Neural Netw., 2010, pp. 1–8. [17] A. Delgado, “The memristor as controller,” in IEEE Nanotechnology Materials and Devices Conference, 2010, pp. 376–379. [18] J. S. A. Doria-Cerezo, L. van der Heijden, “Memristive port-hamiltonian control: Path-dependent damping injection in control of mechanical systems,” in Proc. 4th IFAC Workshop on Lagrangian and Hamiltonian Methods for Non Linear Control, 2012, pp. 167–172. [19] R. Pasumarthy, G. Saha, F. Kazi, and N. Singh, “Energy and power based perspective of memristive controllers,” in Proc. IEEE Conf. Decis. Contr., 2013, pp. 642–647. [20] J. J. Gorman and B. Shapiro, Feedback control of mems to atoms. Springer, 2011. [21] E. Ozalevli and P. E. Hasler, “Tunable highly linear floating-gate cmos resistor using common-mode linearization technique,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 55, pp. 999–1010, 2008. [22] L. O. Chua and S. M. Kang, “Memristive devices and systems,” Proc. IEEE, vol. 64, no. 2, pp. 209–223, 1976. [23] T. A. Wey and W. D. Jemison, “Variable gain amplifier circuit using titanium dioxide memristors,” IET Circuits, Devices & Syst., vol. 5, no. 1, pp. 59–65, 2011. [24] S. Haykin, Communication systems, 5th ed. Wiley Publishing, 2009. [25] S. P. Boyd, Linear matrix inequalities in system and control theory. SIAM, 1994. [26] A. Ben-Tal and A. Nemirovski, Lectures on modern convex optimization: analysis, algorithms, and engineering applications. SIAM, 2001. [27] H. H. Rosenbrook, “The stability of linear time-dependent control systems,” Int. J. Electron. & Contr., vol. 15, pp. 73–80, 1963. [28] V. Solo, “On the stability of slowly time-varying linear systems,” Math. of Contr., Signals & Syst., vol. 7, pp. 331–350, 1994. [29] C. B. Barber, D. P. Dobkin, and H. Huhdanpaa, “The quickhull algorithm for convex hulls,” ACM Trans. on Mathematical Software, vol. 22, no. 4, pp. 469–483, 1996. [Online]. Available: http://www.qhull.org [30] A. Hassibi, J. How, and S. Boyd, “A path-following method for solving bmi problems in control,” in Proc. Amer. Contr. Conf., 1999, pp. 1385– 1389. [31] Y.-Y. Cao, J. Lam, and Y.-X. Sun, “Static output feedback stabilization: an ilmi approach,” Automatica, vol. 34, no. 12, pp. 1641–1645, 1998. [32] K.-C. Goh, M. Safonov, and G. Papavassilopoulos, “A global optimization approach for the bmi problem,” in Proc. IEEE Conf. Decis. Contr, vol. 3, 1994, pp. 2009–2014. [33] K. H. H. Fujioka, “Bounds for the bmi eingenvalue problem - a good lower bound and a cheap upper bound,” Trans. Society of Instr. & Contr. Eng., vol. 33, pp. 616–621, 1997. [34] M. Michihiro Kawanishi and Y. Shibata, “Bmi global optimization using parallel branch and bound method with a novel branching method,” in Proc. Amer. Contr. Conf., 2007, pp. 1664–1669. [35] K. Ogata, Modern Control Engineering, 5th ed. Prentice Hall, 2010. [36] S. Boyd and L. Vandenberghe, Convex optimization. Cambridge university press, 2009. [37] R. T. Rockafellar, Convex analysis. Princeton university press, 1997. [38] A. Hassibi, J. How, and S. Boyd, “Low-authority controller design via convex optimization,” in Proc. IEEE Conf. Decis. Contr., 1998, pp. 1385–1389. [39] I. Barkana, “Simple adaptive control-a stable direct model reference adaptive control methodology-brief survey,” Int. J. Adapt. Contr. Signal Proc., 2013. [40] I. CVX Research, “CVX: Matlab software for disciplined convex programming, version 2.0,” http://cvxr.com/cvx, Aug. 2012. [41] M. K. Mituhiro Fukuda, “Branch-and-cut algorithms for the bilinear matrix inequality eigenvalue problem,” Comput. Optim. Appl, vol. 19, p. 2001, 1999. [42] E. Kreyszig, Introductory functional analysis with applications. Wiley, 1989. Gourav Saha received the B.E. degree from Anna University, Chennai, India, in 2012. Since 2013, he has been working towards his M.S. (by research) degree in the Department of Electrical Engineering at Indian Institute of Technology, Madras. His current research interests include application of memristor and memristive systems, control of MEMS and modelling and control of cyber-physical systems with specific focus on cloud computing. Ramkrishna Pasumarthy obtained his PhD in Systems and Control from the University of Twente, The Netherlands in 2006. He is currently an Assistant Professor at the Department of Electrical Engineering, IIT Madras, India. His research interests lie in the area of modeling and control of physical systems, infinite dimensional systems and control of computing systems, with specific focus on cloud computing systems. Prathamesh Khatavkar received the B.Tech degree from University of Pune, India, in 2012. He is currently pursuing his M.S. (by research) degree in the Department of Electrical Engineering at Indian Institute of Technology, Madras. His current research interests includes Analog/RF IC design.