Tuning of Online Feedback Optimization for setpoint tracking in centrifugal compressors

The controlled system is described by nonlinear dynamics:

where $f:\mathds{R}^{s}\times\mathds{R}^{p}\rightarrow\mathds{R}^{s}$ is continuously differentiable. The inputs $u\in\mathds{R}^{p}$ are constrained by the physics of the system, $-b_{2}\leq u\leq b_{1}$, $b_{i}\in\mathds{R}_{+}^{p}$, $i=1,2$. The outputs $y\in\mathds{R}^{n}$ are described by a continuously differentiable nonlinear mapping $y=g(x,u)$ and $g:\mathds{R}^{s}\times\mathds{R}^{p}\rightarrow\mathds{R}^{n}$. We assume that for a constant $u$, the system (1) reaches a steady state $x_{s}(u)$ such that $f(x_{s},u)=0$. Then we have:

where $h:\mathds{R}^{s}\rightarrow\mathds{R}^{n}$ is a continuously differentiable nonlinear steady state mapping. We further assume that the outputs $y$ in (2) are bounded for any bounded inputs $u$. We want to design an OFO controller so that the outputs (2) track a setpoint $y_{sp}$.

Online Feedback Optimization (OFO) is designed to solve problems of the form (Picallo et al., 2022):

where $\Phi:\mathds{R}^{p}\times\mathds{R}^{n}\rightarrow\mathds{R}$ is a continuously differentiable cost, and $\mathcal{U}=\{u\in\mathds{R}^{p}:Au\leq b\}$, $\mathcal{Y}=\{y\in\mathds{R}^{n}:Cy\leq d\}$, where $A\in\mathds{R}^{q\times p}$, $b\in\mathds{R}^{q}$, $C\in\mathds{R}^{l\times n}$, and $d\in\mathds{R}^{l}$ are constant matrices (Häberle et al., 2020). Online Feedback Optimization (OFO) iteratively updates $u$ in (1) to make $y$ converge to a local optimum of (3), even if $\Phi$ is non-convex (Häberle et al., 2020).

To show the relationship between $\nu$ and $\Delta T$, we analyse the contour plots of the error and the number of oscillations obtained for the constant pressure setpoint, with no setpoint for the torque as functions of the parameters (right column of Fig. 1). The white line indicates where the error and the number of oscillations are equal. We also get that the smallest error is obtained for $\nu=250$ and $\Delta T\leq 0.005$ (bottom right corner in Fig. 1(c)) which corresponds to the largest number of oscillations (bottom right corner in Fig. 1(f)). Conversely, a small $\nu$ and a large $\Delta T$ give no oscillations (top left corner in Fig. 1(f)), but increase the error (top left corner in Fig. 1(c)). The values annotating the intersection line indicate that similar performance can be obtained for $\nu\in[75,175]$ and $\Delta T\in[3,9]$ s, with both error and oscillations approximately 7.2, which suggests that the parameters should be tuned simultaneously.

Figure 1 shows the impact of parameters $\Delta T$ (Fig. 1(a)) and $\nu$ (Fig. 1(d)) on the performance of the OFO controller for a constant setpoint. The oscillatory trajectories in Fig. 1 show the numerical connection between $\nu$ and the derivatives in (8). From (8), we see that:

The parameter $\nu=1$ in Fig. 1(a), so the controller reaches steady state when the optimum is reached and $H^{\top}(u^{k})\nabla\Phi^{\top}(u^{k},y^{k})$ becomes close to zero (vanishing oscillations for $\Delta T\leq 0.05$ in Fig. 1(a)). If $\nu<1$, the derivatives are mitigated, $\nu\|H^{\top}\nabla\Phi^{\top}\|\leq\|H^{\top}\nabla\Phi^{\top}\|$, and the controller increments in (13) are small, leading to sluggish behaviour ($\nu\leq 0.1$ in Fig. 1(d)). Conversely, $\nu>1$ intensifies the effect of the derivatives, $\nu\|H^{\top}\nabla\Phi^{\top}\|\geq\|H^{\top}\nabla\Phi^{\top}\|$, leading to oscillatory behaviour ($\nu\geq 10$ in Fig. 1(d)), in line with the interpretation as a step size in line search algorithms (Bertsekas, 2016, p.31).

To tune OFO, we solve the optimization problem:

The objective (14a) promotes a large sampling time so that the underlying dynamics reaches the steady state within $\Delta T$. The constraint (14b) reinforces the reference tracking in OFO from (11) by restricting the time horizon to $t_{F}$. The constraint (14c) enforces desired properties of the controller without affecting reference tracking. Both constraints are adjusted by choosing $\beta=[\beta_{1},\beta_{2}]$, as suggested by Fig. 1 (right-most column), with $\beta_{1}$ as a threshold for the error and $\beta_{2}$ for the number of oscillations. To solve (14), we use a constant setpoint $p_{sd}^{\text{const}}=0.925$ bar and two trajectories: a truncated sinusoidal signal:

and a step signal

The tuned values from solving (14) were compared with a default OFO controller, with $\Delta T=47.5$ s ensuring timescale separation and corresponding to the settling time of the compressor from (9) (within 5% of the final value (Bagge Carlson et al., 2021)) linearized around the operating points for $\tau=323.6$ Nm.

The system (9) with the controller (7) was simulated using OrdinaryDiffEq.jl (Rackauckas and Nie, 2017). The derivatives $\nabla h$ and $\nabla\Phi$ for solving (5a) were obtained using Zygote v0.6.66 (Innes, 2018). The optimization problem was implemented in an open source programming language, Julia 1.9.0, with Windows (x86_64-w64-mingw32) CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz. The unconstrained optimization problem (5a) was solved using OSQP v0.8.0 (Stellato et al., 2020). To overcome potential non-differentiability (Bachtiar et al., 2016) when solving (14) with respect to the sampling time, we use a derivative-free solver, with a Julia interface NOMAD.jl (Montoison et al., 2020) to NOMAD 4.3.1 (Audet et al., 2022). The parameter $\nu\in[0,10^{3}]$ is bounded for scaling purposes, and $\Delta T\in[5\times 10^{-3},t_{F}/2]$ reflects the physical setup of the compressor defining the smallest and the largest sampling time available. An analysis of the impact of number of iterations $t_{f}/\Delta T$ on the performance of OFO was done by He et al. (2023). We chose $t_{F}=200$ s to reflect typical operating times of centrifugal compressors (Cortinovis et al., 2015).

The results of tuning OFO for the compressor (9) are shown in Fig. 2 with the initial guess for parameters $\Delta T=50$, $\nu=0.1$ in red, the setpoint in dashed black (top row), and the bounds for the torque in black (bottom row). The initial guess for $\nu$ and $\Delta T$ was chosen following the recommendation from Gil et al. (2023), as indicated in Table 1. The values corresponding to the default OFO controller are in solid black lines. To achieve the desired performance, the values of $\beta$ in (14) were applied in decreasing order, starting from $\beta=[150,50]$, as suggested by Fig. 1. Figure 2 shows that for $\beta_{1}=150$ and $\beta_{2}=50$, the constraints on the error and the number of oscillations have little impact on shaping the response ($\Delta T=99$ in Table 1), and lead to a sluggish controller ($\beta_{2}=150$ in Fig. 2(d)). Putting higher priority on the oscillations with $\beta_{2}\leq 25$, and decreasing $\beta_{1}$ from 37.5 to nine allows achieving better tracking with damped oscillations, at the expense of decreased sampling time ($\Delta T\leq 12$ in Table 1). The oscillatory behaviour of the control signal in Fig. 2(d) confirms the importance of preserving the bounds on control signal in (7). Restricting the number of oscillations, $\beta_{2}=12$, yield a less aggressive response (Fig. 2(d)), at the expense of increased error.

To validate the tuning results, we track two trajectories for $p_{sd}$: a sinusoidal trajectory and a step trajectory (dashed black in Fig. 3) reflecting possible setpoints for centrifugal compressors. The three sets of parameters for validation, Set 1, Set 2, and Set 3 from Table 1, correspond to the smallest error obtained for constant setpoint, the smallest number of oscillations with the error below 10 for the step setpoint, and the largest $\Delta T$ and decreased error for sinusoidal tuning trajectory, respectively. The manual tuning for constant setpoint was also chosen because it gave the largest decrease in the error and zero oscillations compared to manual tuning with the step or sinusoidal setpoints. The initial condition for validation was set to $p_{s}(0)=0.91745$ bar, $p_{d}(0)=2$ bar, $m(0)=80$ kg s^-1, and $\omega(0)=700.5$ rad s^-1.

The results of validation are shown in Fig. 3 and Table 2 and confirm that solving (14) to find optimal parameters allows shaping the response of the system. Comparing the impact of the timestep $\Delta T$ on the performance, we see that a large value allows avoiding oscillations (Set 3 and Manual in Fig. 3(a)). Thus, the parameters obtained from (14) for the sinusoidal setpoint (Set 3) are comparable with the manual tuning preserving timescale separation. At the same time, a large $\Delta T$ leads to a sluggish controller and introduces delays (Fig. 3(c)). The performance for the parameters in Set 1 and Set 2 is further defined by the value of $\nu$. The oscillatory behaviour is due to including the setpoint in the calculations of the gradients in (1). From (11), we have $\partial^{2}\Phi/\partial y\partial y_{sp}=-0.02y_{sp}$ and thus a large change in the setpoint $y_{sp}$ corresponds to a large change in the gradient $\partial\Phi/\partial y$ in (8). In Set 1, combining the large change in the gradient due to the setpoint with a large $\nu$ in (4) leads to the controller changing significantly within a small $\Delta T$ and introduces oscillations. At the same time, the parameters from Set 1 allow following the sinusoidal trajectory with the smallest error (first row in Table 2) because there are no abrupt changes in the setpoint trajectory.

The results confirm a nonlinear relationship between the sampling time $\Delta T$ and the step size $\nu$ in (5), suggested by Belgioioso et al. (2022). Instead of iterative tuning by adjusting directly $\Delta T$ and $\nu$, which do not have explicit interpretation in terms of the responses of the system, solving the optimization problem (14) enables adjusting the thresholds $\beta_{1}$ and $\beta_{2}$ directly related to the error and the number of oscillations. As a possible choice of $\beta_{1}$, we can take the error (12) obtained for the initial steady state and a chosen tuning trajectory, $\beta_{1}=\gamma_{1}\int\limits_{0}^{t_{F}}(p_{s}(0)-p_{sd}(\xi))^{2}\mathrm{d}\xi.$ The initial value for $\beta_{2}$ can be chosen with respect to the tuning horizon $t_{F}$, for instance $t_{F}/2$, allowing one zero crossing per two time units. Then both $\beta_{1}$ and $\beta_{2}$ can be decreased until the desired performance is reached.