Encrypted extremum seeking for
privacy-preserving PID tuning as-a-Service

ES is used to determine or maintain an extremum value $J(\boldsymbol{\theta}^{\ast})$ of an objective function $J(\boldsymbol{\theta}):\mathbb{R}^{p}\rightarrow\mathbb{R}$ with parameters $\boldsymbol{\theta}\in\mathbb{R}^{p}$. This is achieved by perturbing $\boldsymbol{\theta}$ at iteration step $k\in\mathbb{N}$ according to

$$\boldsymbol{\theta}_{\boldsymbol{d}}(k):=\boldsymbol{\theta}(k)+\boldsymbol{d}_{\boldsymbol{\theta}}(k)$$

(1)

with $\boldsymbol{d}_{\boldsymbol{\theta}}(k)$ being a user-defined perturbation. Clearly, this results in an “exploration” of $J(\boldsymbol{\theta}_{\boldsymbol{d}}(k))$ and, hence, allows approximating the gradient $\nabla_{\boldsymbol{\theta}}J(\boldsymbol{\theta}(k))$ (or even higher order derivatives) in the vicinity of $\boldsymbol{\theta}(k)$. Based on this information, the parameters are updated by, e.g., a gradient descent step of the form

$$\boldsymbol{\theta}(k+1):=\boldsymbol{\theta}(k)+\Delta\boldsymbol{\theta}(k),$$

(2)

where $\Delta\boldsymbol{\theta}(k):=-\alpha\nabla_{\boldsymbol{\theta}}J(\boldsymbol{\theta}(k))$ and $\alpha\in\mathbb{R}_{+}$ denotes the step width. By repeating these steps for suitable $\alpha$, the algorithm approaches a local minimum $J(\boldsymbol{\theta}^{\ast})$.

Controller tuning. In classical setups, $J(\boldsymbol{\theta})$ is often a cost function that characterizes an optimal operating point, e.g., for power tracking in photovoltaic systems or for maximizing the output of bioreactors [12]. However, ES is also useful for PID tuning [13]. There, the integrated square error

$$J(\boldsymbol{\theta}):=\int e^{2}(t,\boldsymbol{\theta})\mathrm{d}t$$

(3)

is minimized over a user-defined horizon, where $e(t,\boldsymbol{\theta}):=r(t)-y(t,\boldsymbol{\theta})$ is the control error, i.e., the difference between the reference $r(t)$ and the plant output $y(t,\boldsymbol{\theta})$ at time $t\in\mathbb{R}$, and where $\boldsymbol{\theta}=(K_{p}\,\,\,\,K_{i}\,\,\,\,K_{d})^{\top}$ reflects the tuning parameters of the PID

$$C(s,\boldsymbol{\theta}):=K_{p}+K_{i}\frac{1}{s}+K_{d}s.$$

(4)

Remarkably, the extremum seeker has only access to closed-loop information in terms of (3). In fact, additional knowledge about the plant itself is not assumed.

Now, an important difference between classical ES control and ES-based controller tuning is that perturbations $\boldsymbol{\theta}_{\boldsymbol{d}}(k)$ of the controller parameters will not (or only insignificantly) excite the closed-loop in steady state. Hence, additional reference steps of the form

$$r(t)=\begin{cases}\hat{r}\text{ if }t\geq t_{0}\\ 0\text{ else}\end{cases}$$

(5)

are considered in [13] throughout the tuning process. However, other reference signals are possible. Yet, it should be noted that $\nabla_{\boldsymbol{\theta}}J(\boldsymbol{\theta})$ depends on $r(t)$ by means of (3).

Tuning setup. In principle, we have now collected all crucial components for ES-based PID tuning. However, it remains to detail the choice of the perturbations $\boldsymbol{d}_{\boldsymbol{\theta}}(k)$ and the approximation of the gradient $\nabla_{\boldsymbol{\theta}}J(\boldsymbol{\theta}(k))$. The $i\in\{1,\dots,p\}$ perturbation signals $\left(\boldsymbol{d}_{\boldsymbol{\theta}}\right)_{i}$ are in general periodic functions with an amplitude of $\gamma_{i}\in\mathbb{R}$, zero mean, and non-zero squared average (power), e.g., cosine or square waves. With $\boldsymbol{d}_{\boldsymbol{\theta}}$ at hand, the approximation of $\nabla_{\boldsymbol{\theta}}J(\boldsymbol{\theta}(k))$ works as follows. First, the controller parameters are updated by $\boldsymbol{\theta}_{\boldsymbol{d}}(k)$ and $r(t)$ is then used to excite the closed-loop. This allows to measure $e(t,\boldsymbol{\theta})$ over some time frame and to obtain $J(\boldsymbol{\theta}_{\boldsymbol{d}}(k))$ by means of (3). From here on, depending on the choices for $\boldsymbol{d}_{\boldsymbol{\theta}}$, different approaches are possible. For instance, if $\boldsymbol{d}_{\boldsymbol{\theta}}$ are cosine wave perturbations, $\nabla_{\boldsymbol{\theta}}J(\boldsymbol{\theta}(k))$ can be estimated by a high-pass filter and a demodulation signal. This, more classic approach, is used in [13]. Square wave perturbations, on the other hand, enable numerical differentiation [14]. Finally, $\alpha$ and $\gamma_{i}$ control the convergence speed and exploration, respectively. Thus, they must be chosen with care such that the closed-loop does not become unstable during (1) or (2) and the convergence speed is acceptable. For different closed-loop dynamics, these often have to be adapted manually. A convergent controller tuning will, however, improve the controller’s performance by $J(\boldsymbol{\theta})-J(\boldsymbol{\theta}^{\ast})$.

In a nutshell, HE allows carrying out (simple) mathematical operations on encrypted data (see [15] for an introduction). Over the previous decade, various schemes have been proposed with different functionalities (see, e.g, [1, 2, 3, 4]). Commonly supported homomorphisms are encrypted additions and (or) multiplications. More precisely, let $\mathtt{x}_{1}$ and $\mathtt{x}_{2}$ be two arbitrary numbers in the cryptosystem’s plaintext space and let us denote the encryption by “$\mathtt{Enc}$” and the decryption by “$\mathtt{Dec}$”. Then, the homomorphisms “$\otimes$” and “$\oplus$” allow computing

	$\displaystyle\mathtt{x}_{1}\,\mathtt{x}_{2}$	$\displaystyle=\mathtt{Dec}\left(\mathtt{Enc}(\mathtt{x}_{1})\otimes\mathtt{Enc}(\mathtt{x}_{2})\right),$		(6a)
	$\displaystyle\mathtt{x}_{1}+\mathtt{x}_{2}$	$\displaystyle=\mathtt{Dec}\left(\mathtt{Enc}(\mathtt{x}_{1})\oplus\mathtt{Enc}(\mathtt{x}_{2})\right).$		(6b)

Some further operations can easily be derived from these fundamental ones. For instance, we will use “$\ominus$” for encrypted subtractions and “$\boxtimes$” for encrypted multiplications with one public factor $\mathtt{x}_{1}$ or $\mathtt{x}_{2}$. Obviously, with these homomorphisms at hand, polynomial expressions can be evaluated easily. However, non-polynomial expressions are computationally expensive and typically require “tricks”.

Plaintext space and mapping. The plaintext space of HE cryptosystems is usually some finite set with (for example) the canonical representatives

$$\mathbb{Z}_{q}:=\left\{0,\dots,q-1\right\},$$

where $q\in\mathbb{N}_{>1}$ is the ciphertext modulus, i.e., the number of elements in the set $\mathbb{Z}_{q}$. Hence, before one can use HE, all quantities must be mapped to or encoded in $\mathbb{Z}_{q}$. An easy way to do this is as follows. For some scalar $x\in\mathbb{R}$, we compute

$\displaystyle\mathtt{x}:=\lfloor cx\rceil\,\mathrm{mod}\,q,$

(7)

where $c\in\mathbb{R}_{\geq 1}$ is a scaling factor. In order to reconstruct $x$, one can use

$$x\approx\mu(\mathtt{x})/c\;\text{ with }\;\mu(\mathtt{x}):=\begin{cases}\mathtt{x}-q&\text{if }\mathtt{x}\geq q/2\\ \mathtt{x}&\text{else},\end{cases}$$

(8)

where $|x-\mu(\mathtt{x})/c|\leq 1/(2c)$ if $\lfloor cx\rceil\in\mathbb{Z}_{q}$. For multidimensional quantities, i.e., vectors or matrices, (7) and (8) can be used component-wise. Now, it is possible to investigate ciphertext additions and multiplications by means of plaintext additions and multiplications over $\mathbb{Z}_{q}$ due to (6). Then, note that, e.g., $\mu\left(\mathtt{x}^{2}+\mathtt{y}\right)/c^{2}$ is meaningful for $\mathtt{x}$ as in (7) only if $\mathtt{y}:=\lfloor c^{2}y\rceil\,\mathrm{mod}\,q$. In other words, summands should share the same scaling and multiplications increase the scaling. Due to the latter, $q$ and $c$ must be chosen carefully such that overflows, i.e., $\left(\mathtt{x}^{2}+\mathtt{y}\right)\notin\mathbb{Z}_{q}$, are avoided while maintaining sufficient accuracy. Furthermore, for arbitrarily many multiplications with $c>1$ and finite $q$, the plaintext will (almost certainly) overflow $\mathbb{Z}_{q}$, which, e.g., poses a problem for the operating time of a controller [8].

LWE-based cryptosystems. Cryptosystems providing (6) can, e.g., be constructed based on learning with errors (LWE) [16]. These cryptosystems combined with a regular execution of a technique called “bootstrapping” (see [2]) results in a fully homomorphic cryptosystem and resolves the aforementioned problem regarding operating time. Fully HE (theoretically) allows computing any function in terms of boolean or arithmetic circuits. However, bootstrapping is computationally costly (many seconds to minutes up to now in the arithmetic case) such that fully HE is not yet useful for our application. Opposed to fully HE, a leveled HE scheme waives on bootstrapping [but also provides (6)], which results in performance advantages. Consequently, only a finite number of operations¹¹1Additionally to the scale factor growth during multiplications, LWE-based cryptosystems inject a small error in the ciphertexts which grows during operations such that additions and multiplications can only be evaluated finitely often. can be supported, i.e., for a predefined $c$, one must a priori select a large enough $q$ for the arithmetic circuit to be encrypted. Finally, we point out that among the leveled HE schemes, [4] stands out because multiple optimizations such as constructions over polynomial rings with tailored plaintext encodings, residue number systems, or ciphertext rescaling make it somewhat practical. A detailed explanation of the cryptosystem or these methods is, however, out of this scope of this paper.

Several challenges have to be overcome in order to realize encrypted ES-based PID tuning as-a-Service. First (and probably most important), we need to carefully select the parameters of the extremum seeker a priori such that it converges for different unknown closed-loop dynamics without being too conservative. Second, for the estimation of $\nabla_{\boldsymbol{\theta}}J(\boldsymbol{\theta}(k))$ in [13] a high-pass filter is necessary which reuses old output values and would require bootstrapping. Thus, a different technique must be used. Similarly, the integrator for $\boldsymbol{\theta}(k+1)$ can not have unlimited operating time. As specified below, solving these issues will require tailored modifications of the scheme from Section II.A.

Before we move on, let us briefly note that PID controllers of the form (4) are not practical in many industrial applications because high-frequency measurement noise and step-like changes in $r(t)$ can cause very large control inputs due to the derivative part. The standard solution for this issue is to introduce a first order derivative filter with the time constant $T_{f}$. These (more general) PID controllers of the form

$\displaystyle C(s,\boldsymbol{\theta}):=K_{p}+K_{i}\frac{1}{s}+K_{d}\frac{s}{sT_{f}+1}$

with $\boldsymbol{\theta}=\begin{pmatrix}K_{p}\,\,\,\,K_{i}\,\,\,\,K_{d}\,\,\,\,T_{f}\end{pmatrix}^{\top}$ will therefore be considered in the following. Due to the general framework provided by ES, this does not entail changes.

Parameter updates. In order to achieve a convergent ES scheme for different a priori unknown plants, a simple building block is as follows. Instead of absolute updates as in (1) and (2), we make use of

	$\displaystyle\boldsymbol{\theta}_{d\pm}(k)$	$\displaystyle:=\boldsymbol{\theta}(k)\circ\left(\boldsymbol{1}\pm\boldsymbol{d}_{\boldsymbol{\theta}}(k)\right),$		(9a)
	$\displaystyle\boldsymbol{\theta}(k+1)$	$\displaystyle:=\boldsymbol{\theta}(k)\circ\left(\boldsymbol{1}+\Delta\boldsymbol{\theta}(k)\right),$		(9b)

where $\boldsymbol{d}_{\boldsymbol{\theta}}$ and $\Delta\boldsymbol{\theta}$ are relative updates opposed to before. This automatically takes the current size of $\left(\boldsymbol{\theta}(k)\right)_{i}$ into account and results in “sensible” updates. Consider, for example, $K_{p}(k)=100$ and $T_{f}(k)=10^{-3}$. Then, $K_{p}(k+1)=K_{p}(k)(1-0.05)$ and $T_{f}(k+1)=T_{f}(k)(1-0.05)$, i.e., a change by $5\%$, are probably useful updates, whereas $K_{p}(k+1)=K_{p}-0.05$ will likely have a negligible effect and $T_{f}(k+1)=T_{f}(k)-0.05$ results in positive feedback, which may destabilize the closed-loop. Furthermore, we move these updates to the controller. Thus, our tuning scheme will provide $\Delta\boldsymbol{\theta}(k)$ and $\boldsymbol{d}_{\boldsymbol{\theta}}$.

Gradient estimation. Next, we focus on the estimation of $\nabla_{\boldsymbol{\theta}}J(\boldsymbol{\theta})$. Here, instead of using a high-pass filter in combination with a demodulation, we focus on square wave perturbation signals for $\boldsymbol{\theta}_{d\pm}$. In combination with the parameter updates, this allows for arbitrarily many tuning steps and circumvents the necessity for bootstrapping. The result of square wave perturbations is a gradient approximation based on finite differences. However, $2p$ cost function measurements are then necessary in order to approximate $\nabla_{\boldsymbol{\theta}}J(\boldsymbol{\theta})$ which is inefficient. Instead, simultaneous perturbation (see [17]) approximates the (relative) gradient stochastically based on

$$\left(\nabla_{\boldsymbol{\theta}}\tilde{J}(\boldsymbol{\theta})\right)_{i}:=\frac{J\left(\boldsymbol{\theta}_{d+}\right)-J\left(\boldsymbol{\theta}_{d-}\right)}{2\left(\boldsymbol{d}_{\boldsymbol{\theta}}\right)_{i}}$$

(10)

with only two cost function measurements. Here, $\boldsymbol{d}_{\boldsymbol{\theta}}:=\gamma\boldsymbol{h}$ consists of an amplitude $\gamma$ and a random perturbation vector $\boldsymbol{h}$. The entries of $\boldsymbol{h}$ are mutually independent zero-mean random variables, e.g., symmetrically Bernoulli distributed around $0$. Now, using (10) results in various benefits. First, the (expected) performance is higher. More precisely, compared to other approximation techniques, simultaneous perturbations requires the least total amount of cost function evaluations until convergence [14]. Moreover, since $\boldsymbol{\theta}_{d\pm}$ is a square wave signal, faster convergence compared to other perturbation wave forms [18] can be expected. Second, in comparison to [13], the high-pass and additional perturbation signal parameters are now obsolete. Third, it is well-suited for encryption because there are $2^{p}$ possible values for $\boldsymbol{d}_{\boldsymbol{\theta}}$ which allows precomputing the corresponding multiplicative inverses a priori.

Normalization. In the current form, the magnitude of $\nabla_{\boldsymbol{\theta}}\tilde{J}(\boldsymbol{\theta})$ may be very different for different closed-loops, which may also result in quite different $\Delta\boldsymbol{\theta}(k):=\alpha\nabla_{\boldsymbol{\theta}}\tilde{J}(\boldsymbol{\theta})$ and destabilizing or conservative updates. In this context, it would be useful to replace $\nabla_{\boldsymbol{\theta}}\tilde{J}(\boldsymbol{\theta})$ with the normalization $\nabla_{\boldsymbol{\theta}}\tilde{J}(\boldsymbol{\theta})/\|\nabla_{\boldsymbol{\theta}}\tilde{J}(\boldsymbol{\theta})\|$. Then, $\Delta\boldsymbol{\theta}(k)$ could be tuned solely by means of $\alpha$. The problem is, however, that $\|\nabla_{\boldsymbol{\theta}}\tilde{J}(\boldsymbol{\theta})\|$ depends on $\boldsymbol{\theta}(k)$. Consequently, a costly and/or unreliable encrypted inversion would be necessary. This is why we proceed differently here. Our goal is to obtain similar magnitudes for $J(\boldsymbol{\theta}_{d\pm})$ regardless of the closed-loop response. Then, by means of (10) with a predefined $\boldsymbol{d}_{\boldsymbol{\theta}}$, also $\nabla_{\boldsymbol{\theta}}\tilde{J}(\boldsymbol{\theta})$ will be similar in magnitude throughout the tuning. To this end, we replace (3) in (10) by a Newton-Cotes formula

$$\tilde{J}(\boldsymbol{\theta}):=\frac{1}{N}\sum_{n=0}^{N-1}w_{n}\tilde{e}(n,\boldsymbol{\theta})^{2},$$

(11)

where $n\in\mathbb{N}$ denotes the $n$-th discrete time sample of $t$, $\tilde{e}(n,\boldsymbol{\theta}):=1-y(n,\boldsymbol{\theta})/\hat{r}$ with $\hat{r}\neq 0$, $w_{n}$ denote the integration weights, and $N$ refers to the number of samples. Note that the proposed changes amount to scaling a discrete time version of (3) by $1/(T_{\mathrm{int}}\,\hat{r}^{2})$, where $T_{\mathrm{int}}$ is the integration time. Thus, optimal parameter values $\boldsymbol{\theta}^{\ast}$ are not affected. Now, the normalized control error $\tilde{e}(n,\boldsymbol{\theta})$ contributes to a decoupling of (11) from $\hat{r}$. Next, in order to capture the full step response, we propose $N\approx T_{s}/\Delta t$, where $\Delta t$ is the time step and $T_{s}$ is the settling time of the closed-loop. Then, if $N$ is large enough in (11), the magnitude of $\tilde{J}(\boldsymbol{\theta})$ becomes insensitive to variations in $T_{s}$ and $\Delta t$, which naturally happen for different closed-loops.

Survivorship bias. Obviously, the presented tuning approach is somewhat naive because closed-loop stability is not guaranteed during the tuning process as in other model-free tuning methods (see, e.g., [13, 19]). However, as a first step towards privacy-preserving controller tuning, our adaptions of ES work surprisingly well, as we will see in Section IV.

With these preparations and by means of the homomorphisms from Section II.B, an encrypted evaluation of our tailored extremum seeker, in which costly computations are deliberately avoided, can be realized as follows. We first note that, once $\alpha$, $\gamma$, and $c$ are fixed, all $2^{p}$ encrypted perturbations $\mathtt{Enc}\left(\lfloor c\boldsymbol{d}_{\boldsymbol{\theta}}\rceil\,\mathrm{mod}\,q\right)$ and the encrypted multiplicative inverses $\mathtt{Enc}\left(\lfloor(c\alpha)/(2\left(\boldsymbol{d}_{\boldsymbol{\theta}}\right)_{i}\rceil\,\mathrm{mod}\,q\right)$ can be precomputed, where encrypted quantities are scaled in line with (7). Similarly, we precompute $\mathtt{Enc}\left(\lfloor c/\hat{r}\rceil\,\mathrm{mod}\,q\right)$, $\mathtt{Enc}\left(\lfloor c1\rceil\,\mathrm{mod}\,q\right)$, and $\mathtt{Enc}\left(0\right)$, where the latter serves as an initialization for the integration. During operation, $\mathtt{Enc}\left(\lfloor c\boldsymbol{d}_{\boldsymbol{\theta}}\rceil\,\mathrm{mod}\,q\right)$ is randomly selected from the preencrypted values and used as a perturbation signal. Then, after exciting the closed-loop, the cloud obtains $\mathtt{Enc}\left(\mathtt{y}(n)\right)$ and computes $\mathtt{Enc}\left(\tilde{\mathtt{E}}(n)\right)$ as well as $\mathtt{Enc}\left(\tilde{\mathtt{J}}_{j}\right)$ with a step-wise evaluation of (11) for each $j\in\{1,2\}$. For the relative parameter update $\mathtt{Enc}\left(\Delta\text{\straighttheta}\right)$, $\mathtt{Enc}\left(\lfloor(c\alpha)/\left(2\left(\boldsymbol{d}_{\theta}\right)_{i}\right)\rceil\,\mathrm{mod}\,q\right)$ is selected based on the selected perturbation signal. This implementation prohibits the cloud from learning input output data or intermediate results and allows for an arbitrary amount of tuning steps.

Our encrypted scheme (with some additional details) is depicted in Figure 1. There, ciphertext rescaling is applied after every homomorphic (and plaintext) multiplication. We briefly mentioned this technique before in Section II.B. Without losing ourselves in technicalities, with ciphertext rescaling one can reduce the scaling factor of an encrypted plaintext (otherwise it would be $c^{6}$ in $\Delta\text{\straighttheta}$ and not $c$) at the cost of also reducing the modulus (see [20, 4] for details). Thus, the overflow problematic is not resolved. However, apart from some advantages, it is also convenient to use and often automatically performed in popular libraries. An implementation according to Figure 1 can, e.g., be realized with the cryptosystems from [20, 4].

Finally, the remaining parameters are specified as follows. We assume that the number of tuning steps $k_{\max}$, $\hat{r}$, and $\tilde{T}_{s}$ are user-defined, where $\tilde{T}_{s}$ is an estimation of $T_{s}$. We further note that, $N$ as well as $\Delta t$ are public, since an attacker can always count the number of $\mathtt{Enc}(\mathtt{y}(n))$ and measure the time between $\mathtt{Enc}(\mathtt{y}(n+1))$ and $\mathtt{Enc}(\mathtt{y}(n))$. Consequently, also $\tilde{T}_{s}=N\Delta t$ and $y(N)\approx y(N+l)$ for some $l\in\mathbb{N}\setminus\{0\}$ are public.²²2In fact, this enables an encrypted estimation of $\tilde{T}_{s}$ based on comparisons of $\mathtt{Enc}(\mathtt{y}(n))$.