# New Measures for Shaping Trajectories
in Dynamic Optimization

Joshua L. Pulsipher  Benjamin R. Davidson  and Victor M. Zavala  Department of
Chemical and Biological Engineering,
University of Wisconsin-Madison, Madison, WI 53706 USA
(e-mail: {pulsipher,victor.zavala}@wisc.edu).

###### Abstract

We propose a new class of measures for shaping time-dependent trajectories in
dynamic optimization (DO). The proposed measures are analogous to risk
measures used in stochastic optimization (SO) and are inspired by a recently-
proposed unifying abstraction for infinite-dimensional optimization. Risk
measures are summarizing statistics (e.g., average, variance, quantiles,
worst-case values) that are used to shape the probability density of random
objectives and constraints. We show that this extensive collection of measures
can be applied in DO for computing and manipulating interesting features of
time-dependent trajectories (e.g., excursion costs and quantiles). We also
discuss how to implement these measures in the Julia modeling package
InfiniteOpt.jl.

###### keywords:

Dynamic Optimization, Measures, Trajectories, Infinite-Dimensional
Optimization

††thanks: This work was supported by the U.S. Department of Energy under grant
DE-SC0014114.

##  1 Introduction

We consider the continuous-time dynamic optimization (DO) problem:

$\displaystyle\min_{y(\cdot)\in\mathcal{Y}}$ $\displaystyle
M_{t}f(\dot{y}(t),y(t),t)$ (1) s.t. $\displaystyle g(\dot{y}(t),y(t),t)\leq
0,$ $\displaystyle t\in\mathcal{D}_{t}.$

This representation follows from a recently-proposed abstraction that unifies
infinite-dimensional optimization (InfiniteOpt) problems (e.g.,
dynamic/stochastic/PDE optimization) (Pulsipher et al., In Press 2021). In the
context of DO, the infinite domain $\mathcal{D}_{t}$ is the time domain
$[t_{0},t_{f}]$, $y(t)\in\mathcal{Y}\subseteq\mathbb{R}^{n_{y}}$ are time-
valued decision functions (e.g., state/control variables),
$\dot{y}(t)\in\mathbb{R}^{n_{y}}$ are derivative variables $dy(t)/dt$,
$f(\cdot)$ is an infinite-dimensional cost function, and $g(\cdot)$ is a
vector-valued infinite-dimensional constraint functions
$g_{j}(\cdot),\;j\in\mathcal{J}\subseteq\mathbb{R}^{n_{g}}$.

The measure operator $M_{t}:\mathcal{D}_{t}\mapsto\mathbb{R}$ is the focus of
this work; this operator seeks to scalarize infinite-dimensional cost
functions (summarizing them over the time domain $\mathcal{D}_{t}$) to form a
well-posed objective function. Such measures can also be used to handle
constraints, but here we focus on objectives to simplify the presentation.
Problem (1) is general and captures a wide range of DO problems (e.g., model
predictive control and state/parameter estimation) where $y(\cdot)$ is
comprised of state/control/parameter variables and $g(\cdot)\leq 0$ can
include DAE/path/point constraints.

Classical DO formulations minimize the integral of the cost trajectory (known
as the Bolza objective). In other words, these formulations use the measure
operator:

$M_{t}f(t)=\int_{t\in\mathcal{D}_{t}}f(t)dt$ (2)

where we write $f(t):=f(\dot{y}(t),y(t),t)$ for convenience. Minimizing (2)
has the effect of uniformly shaping the cost trajectory $f(t)$ over the domain
$\mathcal{D}_{t}$. This also amounts to minimizing the total cost (which is
equivalent to the average cost with a normalization factor of
$(t_{f}-t_{0})^{-1}$). However, we can envision more advanced DO formulations
that aim to shape the cost trajectory $f(\cdot)$ with other measures (e.g.,
peak or excursion costs). For example, Risbeck and Rawlings recently proposed
an MPC objective that penalizes the total and peak of the cost (Risbeck and
Rawlings, 2019). Minimizing the peak of a cost-trajectory is often a desirable
feature in systems, as these can be associated with unsafe behavior and
economic penalties.

Under the InfiniteOpt representation, Problem (1) is analogous to the SO
problem:

$\displaystyle\min_{z\in\mathcal{Z},y(\cdot)\in\mathcal{Y}}$ $\displaystyle
M_{\xi}f(z,y(\xi),\xi)$ (3) s.t. $\displaystyle g(z,y(\xi),\xi)\leq 0,$
$\displaystyle\xi\in\mathcal{D}_{\xi}.$

where $z\in\mathcal{Z}\subseteq\mathbb{R}^{n_{z}}$ are here-and-now decision
variables, $\xi$ is a random parameter,
$\mathcal{D}_{\xi}\subseteq\mathbb{R}^{n_{\xi}}$ is the co-domain of its
distribution, and $y(\xi)\in\mathcal{Y}\in\mathbb{R}^{n_{y}}$ are recourse
decision variables. We observe that (1) is a special case of (3) if we define
$n_{z}=0$, $n_{\xi}=1$, $\mathcal{D}_{\xi}=\mathcal{D}_{t}$, $M_{\xi}=M_{t}$,
and set $y(\xi)$ to contain both $\dot{y}(t)$ and $y(t)$.

In SO, the measure $M_{\xi}$ can draw from a wide collection of risk measures
(summarizing statistics) that aim to shape the probability density function of
$f(\xi)$; some examples include the expected value, mean-variance, absolute
deviation, and conditional-value-at-risk (CVaR) (Ruszczyński and Shapiro,
2006). Risk measures enable greater flexibility in shaping the cost function
$f(\xi)$ and are commonly used to penalize extreme events (i.e., high costs).
By leveraging the connection between (1) and (3), we have recently shown that
the expected value and CVaR measures can be used as measures in DO to shape
the time-dependent cost trajectory (Pulsipher et al., In Press 2021).

In this work, we formalize the use of risk measures to DO problems that follow
(1). This new class of dynamic measures $M_{t}$ enable us to shape
trajectories and policies in new ways. Moreover, we highlight how the Julia-
based modeling package InfiniteOpt.jl provides a useful interface to compactly
express such measures.

##  2 Main Results

###  2.1 Time-Valued Density Functions

We beging by establishing time analogues of probability density functions
(pdf) and of cumulative density functions (cdf) typically used in SO (i.e.,
pdfs and cdfs defined over the time domain $\mathcal{D}_{t}$). These provide
key constructs to facilitate the use of risk measures in DO.

In SO, $\xi$ is described by the pdf
$p_{\xi}:\mathcal{D}_{\xi}\mapsto\mathbb{R}_{\geq 0}$, which satisfies
$\int_{\xi\in\mathcal{D}_{\xi}}p_{\xi}(\xi)d\xi=1$. The pdf is used to compute
statistics such as the expectation:

$\mathbb{E}_{\xi}[f(\xi)]=\int_{\xi\in\mathcal{D}_{\xi}}f(\xi)p_{\xi}(\xi)d\xi.$
(4)

This measure $M_{\xi}=\mathbb{E}_{\xi}$ is commonly used in SO; here, we can
see that $p_{\xi}(\cdot)$ acts as a weighting function that places varied
emphasis over the domain $\mathcal{D}_{\xi}$. The expectation summarizes the
cost function in a single scalar value.

Following the analogy of SO and DO, we define a weighting function
$p_{t}:\mathcal{D}_{t}\mapsto\mathbb{R}_{\geq 0}$ (i.e., a time-valued pdf) to
yield the time average (expectation):

$\mathbb{E}_{t}[f(t)]:=\int_{t\in\mathcal{D}_{t}}f(t)p_{t}(t)dt.$ (5)

This provides us flexibility in prioritizing different regimes in the domain
$\mathcal{D}_{t}$. We obtain (2) as a special case by setting $p(t)=1$;
however, we can envision choosing from a wide range of candidate pdfs, such
as:

$p_{t}(t)=\frac{1}{t_{f}-t_{0}}$ (6)

Note that this pdf of the uniform random parameter and places equal emphasis
on different parts of the time domain. This inspires the consideration of
other weighting functions, such as the exponential pdf:

$p_{t}(t)=\gamma\e^{-\gamma t}$ (7)

where $\gamma\in\mathbb{R}_{>0}$ is the decay rate parameter. Using this
weighting function to compute the time average places emphasis on the
beginning of the domain and quickly decays over time. In other words, this can
be seen as a discount factor (Shin and Zavala, 2020).

The cdf $P:\mathcal{D}_{\xi}\mapsto[0,1]$ of a random variable $\xi$ is:

$P(\xi;\hat{\xi}):=\mathbb{P}_{\xi}(\xi\leq\hat{\xi})=\int_{\xi\in\\{\xi\in\mathcal{D}_{\xi}:\xi\leq\hat{\xi}\\}}p_{\xi}(\xi)d\xi$
(8)

This denotes the cumulative probability of finding $\xi$ below the threshold
$\hat{\xi}\in\mathcal{D}_{\xi}$. In the context of the DO, we can write the
cdf of the time trajectory $f(t)$ by using the excursion sets:

$\begin{gathered}\mathcal{D}_{t}^{+}(f(t);\hat{f}):=\\{t\in\mathcal{D}_{t}:f(t)\geq\hat{f}\\}\\\
\mathcal{D}_{t}^{-}(f(t);\hat{f}):=\\{t\in\mathcal{D}_{t}:f(t)\leq\hat{f}\\}\end{gathered}$
(9)

where $\mathcal{D}_{t}^{+}(f(t);\hat{f})\subseteq\mathcal{D}_{t}$ and
$\mathcal{D}_{t}^{-}(f(t);\hat{f})\subseteq\mathcal{D}_{t}$ are the positive
and negative function excursion sets, respectively. We can use the negative
excursion set to establish:

$P(f(t);\hat{f})=\int_{t\in\mathcal{D}_{t}^{-}(f(t);\hat{f})}p_{t}(t)dt.$ (10)

In a DO context, the cdf measures the fraction of time that the trajectory
$f(t)$ is below the threshold $\hat{f}$. We will see that expressing the cdf
in terms of excursion set allows us to interpret measures as mechanisms to
bound trajectories.

###  2.2 Dynamic Measures

Here we illustrate the interpretation of risk measures in a DO context. A
large number of risk measures have been proposed in the SO community and
analyzing them is beyond the scope of this work (we refer the reader to
(Krokhmal et al., 2013) for a review). This section establishes key constructs
and exemplifies the steps needed to interpret risk measures in a time setting.
We discuss properties of these proposed measures in Section 2.3.

####  2.2.1 Expectation

The time expectation measure $\mathbb{E}_{t}$ shown in Equation (5) allows us
to assess a weighted average of our cost trajectory $f(t)$ in accordance with
the weighting function $p_{t}(t)$. In the context of DO, minimizing the
expected cost provides an easily interpretable objective and the key modeling
choice lies in the selection of $p_{t}(t)$. Reasonable candidates for many DO
applications are (6) and (7), as demonstrated in Section 3; however, other
weighting functions are possible (e.g., Gaussian and Gamma). The selection of
the weighting function dictates how much emphasis is placed on different parts
of the time domain. The expectation serves as a core construct in defining the
more sophisticated measures.

###### Proposition 1

The time expectation $\mathbb{E}_{t}$ in (5) is a special case of
$\mathbb{E}_{\xi}$ in (4) if $n_{z}=0$, $n_{\xi}=1$,
$\mathcal{D}_{\xi}=\mathcal{D}_{t}$.

####  2.2.2 Mean-Variance

The mean-variance $\mathbb{E}$-$\mathbb{V}_{\xi}$ is a classical measure used
in portfolio optimization (Leland, 1999). For SO problems (3), this minimizes
the variance (i.e., spread) of the cost outcomes in an attempt to mitigate
high-cost events:

$\mathbb{E}\text{-}\mathbb{V}_{\xi}[f(\xi)]:=\mathbb{E}_{\xi}[f(\xi)]+\lambda\mathbb{V}_{\xi}[f(\xi)]$
(11)

where
$\mathbb{V}_{\xi}=\mathbb{E}_{\xi}[(f(\xi)-\mathbb{E}_{\xi}[f(\xi)])^{2}]$ is
the variance and $\lambda\in\mathbb{R}_{\geq 0}$ is a tradeoff parameter.
Transferring this to a DO setting we obtain the time-valued measure:

$\mathbb{E}\text{-}\mathbb{V}_{t}[f(t)]:=\mathbb{E}_{t}[f(t)]+\lambda\mathbb{V}_{t}[f(t)].$
(12)

Minimizing $\mathbb{E}\text{-}\mathbb{V}_{t}[f(t)]$ in a DO problem provides a
tradeoff problem that seeks to minimize the magnitude of the cost trajectory
(the expectation) and the variability/fluctuations (variance) of the cost
trajectory. A disadvantage of this measure is that it penalizes cost
variability equally for low and high costs. However, this property can be
advantageous for cost functions that seek to enforce smooth control
trajectories.

###### Proposition 2

The measure operator $\mathbb{E}$-$\mathbb{V}_{t}$ from (12) is a special case
of $\mathbb{E}$-$\mathbb{V}_{\xi}$ from (11) under the same conditions of
Proposition 1.

####  2.2.3 Quantile

The quantile $Q_{\xi}(f(\xi);\alpha)$ (also referred to as the value-at-risk)
denotes the threshold value $\hat{f}$ for $f(\xi)$ such that the cumulative
probability of incurring costs below the threshold is at least
$\alpha\in[0,1]$:

$Q_{\xi}(f(\xi);\alpha):=\inf_{\hat{f}\in\mathbb{R}}\left\\{P(f(\xi);\hat{f})\geq\alpha\right\\}.$
(13)

Constraining (13) is equivalent to enforcing a probabilistic constraint
(Sarykalin et al., 2008):

$Q_{\xi}(f(\xi);\alpha)\leq 0\iff\mathbb{P}_{\xi}(f(\xi)\leq 0)\geq\alpha.$
(14)

We can use the cdf (10) in combination with (13) to define the time-valued
quantile:

$Q_{t}(f(t);\alpha):=\inf_{\hat{f}\in\mathbb{R}}\left\\{\int_{t\in\mathcal{D}_{t}^{-}(f(t);\hat{f})}p_{t}(t)dt\geq\alpha\right\\}.$
(15)

Using this measure in (1) minimizes the excursion threshold of the cost
function trajectory such that the fraction of time that exceed it is no more
than $1-\alpha$. Unlike the $\mathbb{E}$-$\mathbb{V}_{t}$ measure, the
quantile measure only penalizes high cost values, making it an attractive
alternative in certain cases. However, the potential disadvantages of this
measure are that it does not strongly discourage high cost peaks in the
positive function excursion $\mathcal{D}_{t}^{+}(f(t);Q_{t}(f(t);\alpha))$ and
it is nonconvex and difficult to compute in general.

###### Proposition 3

The quantile measure $Q_{t}$ (15) is a special case of its analogue $Q_{\xi}$
in (13) under the same conditions of Proposition 1.

####  2.2.4 Conditional-Value-at-Risk

The conditional-value-at-risk (CVaR) measure seeks to address the limitations
of the quantile measure $\mathbb{E}$-$\mathbb{V}_{\xi}$ by penalizing the
expected value of the $1-\alpha$ largest cost values:

$\text{CVaR}_{\xi}(f(\xi);\alpha):=\min_{\hat{f}\in\mathbb{R}}\left\\{\hat{f}+\frac{1}{1-\alpha}\mathbb{E}_{\xi}[f(\xi)-\hat{f}]_{+}\right\\}$
(16)

where
$\mathbb{E}_{\xi}[f(\xi)-\hat{f}]_{+}:=\mathbb{E}_{\xi}[\max(0,f(\xi)-\hat{f})]$
and $\alpha\in[0,1)$. CVaR is also known as the superquantile; under mild
assumptions, CVaR can be represented as:

$\text{CVaR}_{\xi}(f(\xi);\alpha)=\mathbb{E}_{\xi}[f(\xi):f(\xi)\geq
Q_{t}(f(t);\alpha)]$ (17)

since the minimizer $\hat{f}^{*}$ is $Q_{t}(f(t);\alpha)$ (Rockafellar et al.,
2000). With this observation, the time-valued CVaR measure can be expressed
as:

$\text{CVaR}_{t}(f(t);\alpha):=\min_{\hat{f}\in\mathbb{R}}\left\\{\hat{f}+\frac{1}{1-\alpha}\mathbb{E}_{t}[f(t)-\hat{f}]_{+}\right\\}.$
(18)

This provides a convex measure that penalizes the high (peak) costs incurred
in the positive function excursion set
$\mathcal{D}_{t}^{+}(f(t);Q_{t}(f(t);\alpha))$. Note that this penalizes
multiple peak costs (and not just the peak cost, as done in typical DO
formulation). Moreover, one can show that:

$\displaystyle\lim_{\alpha\rightarrow
0}\text{CVaR}_{t}(f(t);\alpha)=\mathbb{E}_{t}[f(t)]$ (19)
$\displaystyle\lim_{\alpha\rightarrow
1}\text{CVaR}_{t}(f(t);\alpha)=\max_{t\in\mathcal{D}_{t}}f(t)$

which both follow from Equation (17). As such, CVaR is highly versatile
measure for use in DO that can capture both average and extreme features of a
time trajectory.

###### Proposition 4

The CVaR measure in Equation (18) is a special case of CVaRξ in Equation (16)
under the same conditions of Proposition 1.

####  2.2.5 Disutility

Disutility risk measures are another prevalent measure class used in SO; this
employs an expectation over a disutility function
$g:\mathbb{R}\mapsto\mathbb{R}$ (typically a convex increasing function) that
penalizes unfavorable values of $f(\xi)$:

$D_{\xi}(f(\xi)):=\mathbb{E}_{\xi}[g(f(\xi))].$ (20)

In an effort to make a translation invariant measure, the measure is often
expressed as:

$\tilde{D}_{\xi}(f(\xi)):=\inf_{\hat{f}\in\mathbb{R}}\mathbb{E}_{\xi}[f(\xi)+g(f(\xi)-\hat{f})].$
(21)

One can show that $\text{CVaR}_{\xi}(f(\xi);\alpha)$ is a special case of
(21). This follows by letting $g(x;\alpha)=(1-\alpha)^{-1}\max(0,x)-x$ where
$x\in\mathbb{R}$. Then by substituting $g(x;\alpha)$ in (21) we obtain (16).
We transfer (21) to DO by using the time-valued expectation $\mathbb{E}_{t}$:

$\tilde{D}_{t}(f(t)):=\inf_{\hat{f}\in\mathbb{R}}\mathbb{E}_{t}[f(t)+g(f(t)-\hat{f})].$
(22)

This measure class provides great flexibility in shaping dynamic trajectories
as there are diverse choices of $g(\cdot)$ (in addition to the flexibility
provided via selecting the weighting function $p_{t}(\cdot)$). A useful survey
on the properties of disutility functions in the context of SO is provided in
(Fulga, 2016).

###### Proposition 5

The time-valued disutility measure $\tilde{D}_{t}(f(t))$ is a special case of
$\tilde{D}_{\xi}(f(\xi))$ under the same conditions of Proposition 1.

###  2.3 Measure Properties

Here we formalize some key mathematical properties of the dynamic measures
presented in Section 2.2. This provides some interesting and useful insights
on the behavior that these measures induce. These properties have been studied
in the SO community and we will show that this rich theory can be readily
applied to DO.

In the context of SO, four main properties are typically considered for risk
measures: convexity, monotonicity, translation invariance, and positive
homogeneity. Moreover, a measure operator is said to be coherent if it
satisfies all these properties (Artzner et al., 1999; Ruszczyński and Shapiro,
2006).

Convexity asserts that a measure operator $M_{\xi}$ satisfy:

$M_{\xi}(\beta f+(1-\beta)h)\leq\beta M_{\xi}(f)+(1-\beta)M_{\xi}(h)$ (23)

for all measurable functions
$f(\xi),h(\xi):\mathcal{D}_{\xi}\mapsto\mathbb{R}$ in the linear function
space $\mathscr{F}$ and all $\beta\in[0,1]$. This property is key for creating
optimization objectives that are well-posed and guarantees that the measure of
a convex cost is also convex.

Under monotonicity, we have that if $f_{1}(\xi)\succeq f_{2}(\xi)$
($f_{1}(\xi)$ dominates $f_{2}(\xi)$), then the measure $M_{\xi}$ satisfies:

$M_{\xi}(f_{1}(\xi))\geq M_{\xi}(f_{2}(\xi)).$ (24)

This ensures that, if a cost function dominates another cost function, then
the measure former will also be greater. than that of the latter. The concept
of dominance (comparing whether a random variable is better than another
random variable) is an interesting and important concept that has not been
explored in DO. In a DO context, dominance of first-order ($f_{1}(t)\succeq
f_{2}(t)$) requires that $P(f_{1}(t)>\hat{f})\geq P(f_{2}(t)>\hat{f})$ for any
threshold value $\hat{f}$. In other words, the fraction of time that the
trajectory $f_{1}(t)$ remains above the threshold is greater or equal than the
fraction of time that the trajectory $f_{2}(t)$ remains above the same
threshold. A monotonic measure is such that, if $f_{1}(t)\succeq f_{2}(t)$,
then $M_{t}(f_{1}(t))\geq M_{t}(f_{2}(t))$ holds. Note that dominance holds
trivially if $f_{1}(t)\geq f_{2}(t)$ for all $t\in\mathcal{D}_{t}$. These
concepts are important because comparisons (benchmarks) of time trajectories
are not as straightforward (Renteria et al., 2018), as the trajectories are
functions (not scalar values) and thus a trajectory might be better in some
parts of the time domain but not in others

A translation invariant measure satisfies:

$M_{\xi}(f(\xi)+a)=M_{\xi}(f(\xi))+a$ (25)

if $a\in\mathbb{R}$ and $f(\xi)\in\mathscr{F}$. In a DO context, this property
ensures that offsetting the cost function will not change the shape of the
optimal cost trajectory.

The positive homogeneity property is given by:

$M_{\xi}(\tau f(\xi))=\tau M_{\xi}(f(\xi))$ (26)

if $\tau>0$ and $f(\xi)\in\mathscr{F}$. In the context of DO, a positive
homogeneous measure provides the property that uniformly scaling the cost by
$\tau$ will not affect the shape of the optimal cost trajectory.

The analysis of the stochastic risk measures featured in Section 2.2 is well-
established in the SO literature and Table 1 provides a summary of these
(Artzner et al., 1999; Ruszczyński and Shapiro, 2006).

$M_{\xi}$ | (23) | (24) | (25) | (26)
---|---|---|---|---
$\mathbb{E}_{\xi}$ | Yes | Yes | Yes | Yes
$\mathbb{E}$-$\mathbb{V}_{\xi}$ | Yes | No | Yes | No
$Q_{\xi}$ | No | Yes | Yes | Yes
CVaRξ | Yes | Yes | Yes | Yes
$\tilde{D}_{\xi}$ | Yes | Yes | Yes | Yes
Table 1: A summary of the properties satisfied by certain measures $M_{\xi}$.
Note that $\tilde{D}_{\xi}$ only satisfies (26) if $g(\cdot)$ is positive
homogeneous.

Since the time-valued measures are special cases of the SO counterparts, they
inherit the properties of Table 1. This illustrates how one can transfer rich
theory from SO (with respect to these measure operators) to a DO context. The
time-valued expectation measure is a coherent risk measure; this might explain
why this has been the classical measure used in DO. It is particularly
important to observe that convexity ensures that the use of this measure
yields a convex objective if the cost function is convex. Interestingly, the
convexity of the objective has key implications for establishing stability
conditions (e.g., closed-loop stability of MPC) (Rawlings et al., 2017). From
this, we observe that other non-convex measures such as
$\mathbb{E}$-$\mathbb{V}_{t}$ and $Q_{t}$ may not yield stability. On the
other hand, CVaRt and $\tilde{D}_{t}$ are convex and thus might inherit
stability properties (this is an interesting topic of future work).

###  2.4 Modeling in InfiniteOpt.jl

The unifying abstraction for infinite-optimization (which facilitated the
connection between Problems (1) and (3)) is implemented in a Julia package
called InfiniteOpt.jl (Pulsipher et al., In Press 2021). This enables us to
intuitively model continuous DO problems (e.g., Problem (1)) following a
simple symbolic syntax (see Code Snippet 1). It is measure-centric and readily
enables to quickly implement new candidate measure operators (such as those
proposed in Section 2.2). Moreover, the unifying abstraction behind
InfiniteOpt.jl facilities the incorporation of random constructs and/or PDE
constraints. By default, these models are solved via direct transcription, but
other solution methodologies can be implemented.

##  3 Case Study

We compare the time-valued measures proposed in Section 2.2 in the context of
optimal control. We adapt the pandemic control problem that seeks to choose an
isolation policy to control the spread of a contagion that minimally impacts
the economic impact (imposed by mandated isolation). We model the spread of
the disease via the SEIR model which defines the populations of individuals
susceptible to infection $y_{s}:\mathcal{D}_{t}\rightarrow[0,1]$, exposed
individuals that are not yet infectious
$y_{e}:\mathcal{D}_{t}\rightarrow[0,1]$, infectious individuals
$y_{i}:\mathcal{D}_{t}\rightarrow[0,1]$, and recovered individuals
$y_{r}:\mathcal{D}_{t}\rightarrow[0,1]$ (considered immune to future
infection). Moreover, these satisfy $y_{s}(t)+y_{e}(t)+y_{i}(t)+y_{r}(t)=1$.
Thus, our state variables are comprised of $y_{s}(t)$, $y_{e}(t)$, $y_{i}(t)$,
and $y_{r}(t)$. Moreover, we exhibit control by imposing an isolation policy
$y_{u}(t)\in[0,\overline{y_{u}}]\subseteq[0,1]$ that entails the separation of
susceptible and exposed individuals ($y_{u}(t)=0$ denotes no separation and
$y_{u}(t)=1$ denotes complete separation). The formulation seeks to minimize
the isolation policy function $y_{u}(t)$ while enforcing that the amount of
infectious individuals $y_{i}(t)$ remains below $\overline{y_{i}}\in(0,1]$:

$\displaystyle\min$ $\displaystyle M_{t}y_{u}(t)$ (27) s.t.
$\displaystyle\frac{dy_{s}(t)}{dt}=\left(y_{u}(t)-1\right)\beta
y_{s}(t)y_{i}(t),\ t\in\mathcal{D}_{t}$
$\displaystyle\frac{dy_{e}(t)}{dt}=\left(1-y_{u}(t)\right)\beta
y_{s}(t)y_{i}(t)-\xi y_{e}(t),\ t\in\mathcal{D}_{t}$
$\displaystyle\frac{dy_{i}(t)}{dt}=\xi y_{e}(t)-\gamma y_{i}(t),\
t\in\mathcal{D}_{t}$ $\displaystyle\frac{dy_{r}(t)}{dt}=\gamma y_{i}(t),\
t\in\mathcal{D}_{t}$ $\displaystyle
y_{s}(0)=s_{0},y_{e}(0)=e_{0},y_{i}(0)=i_{0},y_{r}(0)=r_{0}$ $\displaystyle
y_{i}(t)\leq\overline{y_{i}},\ t\in\mathcal{D}_{t}$ $\displaystyle
y_{u}(t)\in\left[0,\overline{y_{u}}\right],\ t\in\mathcal{D}_{t}$

where $s_{0},e_{0},i_{0},r_{0}\in[0,1]$ are initial conditions and
$\beta,\gamma,\xi\in\mathbb{R}$ are the rates of infection, recovery, and
incubation, respectively, which are specific to the disease in question. Here
we specify initial conditions at $s_{0}=.9999$, $e_{0}=10^{-5}$, and
$i_{0}=r_{0}=0$. The disease parameters are taken to be $\beta=0.727$,
$\gamma=0.303$, and $\xi=0.3$. We choose limits $\overline{y_{u}}=0.8$ and
$\overline{y_{i}}=0.02$. Finally, we set $\mathcal{D}_{t}=[0,200]$. We model
Formulation (27) in InfiniteOpt.jl and use backward finite-difference to
evaluate the derivatives using 101 discretization points. Code Snippet 1
highlights the compact syntax required to model Formulation (27) in
InfiniteOpt.jl under these conditions. For the objective $M_{t}y_{u}(t)$, we
consider the following time-valued measures:
$\int_{t\in\mathcal{D}_{t}}y_{u}(t)dt$, $\mathbb{E}_{t}[y_{u}(t)]$,
$\mathbb{E}$-$\mathbb{V}_{t}(y_{u}(t);\lambda)$, and
CVaR${}_{t}(y_{u}(t);\alpha)$. We also investigate the implications of using
the uniform time-valued pdf in (6) against the exponential pdf of (7).

|
---|---
|
Figure 1: The optimal policy trajectories $y_{u}(t)$. Top Left:
$\mathbb{E}_{t}$ with uniform pdf. Top Right: $\mathbb{E}_{t}$ with
exponential pdf ($\gamma=1$). Bottom Left: $\mathbb{E}$-$\mathbb{V}_{t}$ with
$\lambda=8$ and uniform pdf. Bottom Right: CVaRt with $\alpha=0.9$ and uniform
pdf.

⬇

1using$\mathbb{E}$InfiniteOpt,$\mathbb{E}$Ipopt

2

3#$\mathbb{E}$Set$\mathbb{E}$the$\mathbb{E}$parameters

4$\gamma$,$\mathbb{E}$$\beta$,$\mathbb{E}$$\xi$$\mathbb{E}$=$\mathbb{E}$0.303,$\mathbb{E}$0.727,$\mathbb{E}$0.3

5s0,$\mathbb{E}$e0,$\mathbb{E}$i0,$\mathbb{E}$r0$\mathbb{E}$=$\mathbb{E}$1$\mathbb{E}$-$\mathbb{E}$1e-5,$\mathbb{E}$1e-5,$\mathbb{E}$0,$\mathbb{E}$0

6

7#$\mathbb{E}$Define$\mathbb{E}$the$\mathbb{E}$model

8m$\mathbb{E}$=$\mathbb{E}$InfiniteModel(Ipopt.Optimizer)

9

10#$\mathbb{E}$Define$\mathbb{E}$the$\mathbb{E}$time$\mathbb{E}$parameter

11@infinite_parameter(m,$\mathbb{E}$t$\mathbb{E}$$\in$$\mathbb{E}$[0,$\mathbb{E}$200],$\mathbb{E}$num_supports$\mathbb{E}$=$\mathbb{E}$101)

12

13#$\mathbb{E}$Add$\mathbb{E}$the$\mathbb{E}$variables

14@variable(m,$\mathbb{E}$ys,$\mathbb{E}$Infinite(t))

15@variable(m,$\mathbb{E}$ye,$\mathbb{E}$Infinite(t))

16@variable(m,$\mathbb{E}$yi$\mathbb{E}$$\leq$$\mathbb{E}$0.02,$\mathbb{E}$Infinite(t))

17@variable(m,$\mathbb{E}$yr,$\mathbb{E}$Infinite(t))

18@variable(m,$\mathbb{E}$0$\mathbb{E}$$\leq$$\mathbb{E}$yu$\mathbb{E}$$\leq$$\mathbb{E}$0.8,$\mathbb{E}$Infinite(t))

19

20#$\mathbb{E}$Set$\mathbb{E}$the$\mathbb{E}$time$\mathbb{E}$expectation$\mathbb{E}$objective

21@objective(m,$\mathbb{E}$Min,$\mathbb{E}$$\mathbb{E}$$\mathbb{E}$(yu,$\mathbb{E}$t))

22

23#$\mathbb{E}$Define$\mathbb{E}$the$\mathbb{E}$SEIR$\mathbb{E}$equations

24@constraint(m,$\mathbb{E}$$\partial$(ys,$\mathbb{E}$t)$\mathbb{E}$==$\mathbb{E}$-(1$\mathbb{E}$-$\mathbb{E}$yu)$\mathbb{E}$*$\mathbb{E}$$\beta$$\mathbb{E}$*$\mathbb{E}$ys$\mathbb{E}$*$\mathbb{E}$yi)

25@constraint(m,$\mathbb{E}$$\partial$(ye,$\mathbb{E}$t)$\mathbb{E}$==$\mathbb{E}$(1$\mathbb{E}$-$\mathbb{E}$yu)$\mathbb{E}$*$\mathbb{E}$$\beta$$\mathbb{E}$*$\mathbb{E}$ys$\mathbb{E}$*$\mathbb{E}$yi$\mathbb{E}$-$\mathbb{E}$$\xi$$\mathbb{E}$*$\mathbb{E}$ye)

26@constraint(m,$\mathbb{E}$$\partial$(yi,$\mathbb{E}$t)$\mathbb{E}$==$\mathbb{E}$$\xi$$\mathbb{E}$*$\mathbb{E}$ye$\mathbb{E}$-$\mathbb{E}$$\gamma$$\mathbb{E}$*$\mathbb{E}$yi)

27@constraint(m,$\mathbb{E}$$\partial$(yr,$\mathbb{E}$t)$\mathbb{E}$==$\mathbb{E}$$\gamma$$\mathbb{E}$*$\mathbb{E}$yi)

28@constraint(m,$\mathbb{E}$ys(0)$\mathbb{E}$==$\mathbb{E}$ys0)

29@constraint(m,$\mathbb{E}$ye(0)$\mathbb{E}$==$\mathbb{E}$ye0)

30@constraint(m,$\mathbb{E}$yi(0)$\mathbb{E}$==$\mathbb{E}$yi0

31@constraint(m,$\mathbb{E}$yr(0)$\mathbb{E}$==$\mathbb{E}$yr0)

32

33#$\mathbb{E}$Solve$\mathbb{E}$the$\mathbb{E}$model$\mathbb{E}$and$\mathbb{E}$retrieve$\mathbb{E}$results

34optimize!(m)

35u_opt$\mathbb{E}$=$\mathbb{E}$value(yu)

36ts$\mathbb{E}$=$\mathbb{E}$value(t)

Code Snippet 1: Formulation (27) implemented in InfiniteOpt.jl.

Figure 1 summarizes the results; the expectation $\mathbb{E}_{t}[y_{u}(t)]$
with the uniform pdf shapes $y_{u}(t)$ identically to the classical integral
measure. However, we are able to place increased emphasis on the early time
regime when we use the exponential pdf defined in (7) in combination with
$\mathbb{E}_{t}[y_{u}(t)]$. In comparison to the other expectation measure,
the exponentially weighted counterpart exhibits a policy trajectory that is
significantly reduced at early times while later times lead to increased
isolation requirements. This highlights how the choice of pdf $p_{t}(t)$
enhances the flexibility of our proposed measures in accordance with the
requirements of the problem. The optimal policy we obtain with the mean-
variance measure $\mathbb{E}$-$\mathbb{V}_{t}(y_{u}(t);8)$ demonstrates how
placing increased priority on minimizing the variance of the cost function
induces the trajectory to be increasingly smoothed (i.e., cost fluctuations
are damped). This comes at the trade-off (controlled via specification of
$\lambda$) of increasing the mean isolation policy, but helps to derive a more
consistent policy. For this application, a smoother policy would likely be
preferred since rapid policy changes can be highly disruptive and can lead to
public dissatisfaction. Finally, in contrast to the mean-variance (which
equally penalizes positive and negative cost deviations from the mean), the
CVaR measure only penalizes the high cost deviations that surpass the
threshold determined by the $\alpha$-quantile. In Figure 1 we see that
CVaR${}_{t}(y_{u}(t;0.9)$ flattens the peak isolation policy values observed
with the standard integral/expectation policy. This hedging against high costs
also induces a more substantial response at later times which results in a
larger cumulative cost. Thus, we observe a trade-off (controlled via $\alpha$)
between penalizing cost peaks and minimizing the total cumulative cost.

##  4 Conclusions

In this work, we have shown that risk measures used in SO can be transferred
(along with their mathematical properties) to DO. This enables a new class of
DO formulations that shape time trajectories in interesting and useful ways.
The transfer of insights across the SO and DO disciplines is facilitated by a
unifying infinite-dimensional abstraction. In future work, it will be
interesting to investigate the analogy between SO and DO further in
transferring more amenable measure operators and establishing their properties
in a DO context (e.g., dominance and stability). Moreover, the establishment
of time-valued pdfs provides a foundation from which the utility of
transferring distributionally robust measure functions to DO can be
investigated; which would potentially allow us to consider multiple weighting
functions over the time horizon.

## References

  * Artzner et al. (1999) Artzner, P., Delbaen, F., Eber, J.M., and Heath, D. (1999).  Coherent measures of risk.  _Mathematical finance_ , 9(3), 203–228.
  * Fulga (2016) Fulga, C. (2016).  Portfolio optimization with disutility-based risk measure.  _European Journal of Operational Research_ , 251(2), 541–553.
  * Krokhmal et al. (2013) Krokhmal, P., Zabarankin, M., and Uryasev, S. (2013).  Modeling and optimization of risk.  _Handbook of the fundamentals of financial decision making: Part II_ , 555–600.
  * Leland (1999) Leland, H.E. (1999).  Beyond mean–variance: Performance measurement in a nonsymmetrical world (corrected).  _Financial analysts journal_ , 55(1), 27–36.
  * Pulsipher et al. (In Press 2021) Pulsipher, J.L., Zhang, W., Hongisto, T.J., and Zavala, V.M. (In Press 2021).  A unifying modeling abstraction for infinite-dimensional optimization.  _Computers & Chemical Engineering_.
  * Rawlings et al. (2017) Rawlings, J.B., Mayne, D.Q., and Diehl, M. (2017).  _Model predictive control: theory, computation, and design_ , volume 2.  Nob Hill Publishing Madison, WI.
  * Renteria et al. (2018) Renteria, J.A., Cao, Y., Dowling, A.W., and Zavala, V.M. (2018).  Optimal pid controller tuning using stochastic programming techniques.  _AIChE Journal_ , 64(8), 2997–3010.
  * Risbeck and Rawlings (2019) Risbeck, M.J. and Rawlings, J.B. (2019).  Economic model predictive control for time-varying cost and peak demand charge optimization.  _IEEE Transactions on Automatic Control_ , 65(7), 2957–2968.
  * Rockafellar et al. (2000) Rockafellar, R.T., Uryasev, S., et al. (2000).  Optimization of conditional value-at-risk.  _Journal of risk_ , 2, 21–42.
  * Ruszczyński and Shapiro (2006) Ruszczyński, A. and Shapiro, A. (2006).  Optimization of risk measures.  In _Probabilistic and randomized methods for design under uncertainty_ , 119–157. Springer.
  * Sarykalin et al. (2008) Sarykalin, S., Serraino, G., and Uryasev, S. (2008).  Value-at-risk vs. conditional value-at-risk in risk management and optimization.  In _State-of-the-art decision-making tools in the information-intensive age_ , 270–294. Informs.
  * Shin and Zavala (2020) Shin, S. and Zavala, V.M. (2020).  Diffusing-horizon model predictive control.  _arXiv preprint arXiv:2002.08556_.


# The normal Casimir force for lateral moving planes with isotropic
conductivities

N. Emelianova and N. Khusnutdinov111email<EMAIL_ADDRESS>
CMCC, Universidade Federal do ABC, Avenida dos Estados 5001, CEP 09210-580,
SP, Brazil

(October 13, 2023)

###### Abstract

We consider the two planes with isotropic conductivity that are in relative
lateral motion with velocity $v$ and inter-plane distance $a$. Two models of
conductivity are taken into account – the constant and frequency-dependent
Drude model. The normal (perpendicular to planes) Casimir force is analysed in
detail for two systems – i) two planes with identical conductivity and ii) one
of the planes is a perfect metal. The velocity correction to the Casimir
energy $\Delta_{v}\mathcal{E}\sim v^{2}$ for small velocity for all considered
cases. In the case of the constant conductivity $\eta$, the energy correction
is
$\Delta_{v}\mathcal{E}\sim\frac{\eta}{a^{3}}\left(\frac{v}{\eta}\right)^{2}$for
$v\ll\eta\ll 1$.

##  1 Introduction

The Casimir effect [1] was first considered for perfect conductive plates, and
nowadays, it is extended to many non-ideal and new materials [2, 3]. The
relative motions of bodies give an additional contribution to the Casimir
force between bodies (see the recent review on the dynamic Casimir effect [4]
and Refs. [2, 5, 6]). The relative motions are lateral (parallel to the
planes), perpendicular to the planes, or, in general, their combinations. The
Casimir effect for perpendicularly and uniformly moving slabs has been
considered firstly in Refs. [7, 8] for electromagnetic and massless scalar
fields. It is a direct consequence of the Quantum Field Theory with the moving
boundaries [9]. In the non-relativistic case, the velocity correction to the
Casimir pressure is quadratic $\sim v^{2}$ for both fields but with opposite
signs. It is positive for the massless scalar field and it is negative for the
electromagnetic case.

The lateral relative motion of the planes gives two different Casimir
pressures in the perpendicular directions. One of them is normal to the planes
as it is for perpendicular motion, and the second is along planes; it is
called a quantum, non-contact, or Casimir friction. The normal force was
considered in Ref. [10] for layers in stratified dielectric media with
magneto-electric and non-reciprocal coupling. For a system of three layers,
the force may be attractive or repulsive depending on the velocity directions
of extreme layers. In the non-relativistic case, the force becomes repulsive
if extreme layers have the same velocity directions with respect to the middle
layer and it is attractive for opposite directions of velocities. The velocity
correction to the Casimir energy has the same order $\sim v^{2}$. For
relativistic velocities, it may be attractive as well as repulsive.

Quantum friction is more difficult for analysis problems and it is currently
under discussion these days [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
23, 24], up to complete negation of quantum friction existence [16]. Two
dielectric planes at different temperatures with lateral relative motion have
been considered in Ref. [11, 12]. The quantum friction force was calculated in
the framework of the Rytov fluctuation theory [25]. It was shown that it is
proportional to the first power of velocity $v$, but may get different signs –
it may slow down or accelerate planes. Mkrtchian in Ref. [13] considered two
conductive planes with relative lateral motion and calculated the force and
viscosity of vacuum for two models of the plane’s impedance. The dependence of
the force on the inter-plane distance crucially depends on the model selected
but in any case the velocity correction to the usual Casimir force $\sim v$
which is the usual dependence for Casimir friction. The quantum friction for
planes with temporally dispersive conductivity was analysed in Ref. [14]. For
the case of constant conductivity the quantum friction $\sim v^{3}$ at small
velocity and $\sim v^{-1}\ln v$ for high velocity. Volokitin and Persson [15]
developed a general theory for quantum friction in the framework of
fluctuation electrodynamics. The absence of quantum friction as a whole was
claimed in Ref. [16], but the normal force exists (the discussion see in Ref.
[17]). The quantum friction in the framework of scattering theory was
considered in Ref. [18, 19]. It was shown that the quantum friction threshold
exists – it is zero if the relative velocity is smaller than the speed of
light in the materials of slabs. The origin of quantum friction was connected
to the quantum Cherenkov radiation: the super-luminally moving object
spontaneously emits photons. It is closely connected with super-radiance: the
rotating body amplifies incident waves [26]. In Refs. [20, 21], quantum
friction was calculated for two graphene sheets in the framework of the
effective action approach. It was noted the threshold for velocity – the
friction is zero if the relative velocity is smaller than the velocity of
Fermi. It is correlated with [18, 19]: the Dirac electron in graphene is
described by the Dirac equation with velocity Fermi instead of the velocity of
light. The threshold was confirmed by different calculations in Ref. [24]. The
quantum friction in the framework of the nonperturbative approach was
considered in Ref. [22]. The friction force is associated with electromagnetic
instability – the kinetic energy of the relative motion is transformed into an
exponentially growing coherent radiation.

In Ref. [24], the general approach was developed for two conductive planes
with relative lateral motion. In the framework of this approach, the normal
and tangential forces may be calculated. The normal force is reduced to the
usual Casimir force for planes with tensorial conductivities [27] where the
specific form for the tensor of moving plane was used. Quantum friction may
appear as an imaginary part of the energy calculated for some complex
frequencies, which are noted in earlier papers Refs. [28, 29, 30].

In this paper, we consider the normal force for lateral moving planes with
isotropic conductivities. As noted in Refs. [31, 32], the Ohm law for the
moving plane should be considered carefully. The Lorentz transformation even
for a scalar constant conductivity is not simple because it is a coefficient
between $3$–vectors of the electric current and an electric field. The
transformation law of conductivity tensor was discussed in Refs. [32] by using
a linear response tensor [33]. In the case of graphene, the role of linear
response tensor plays the polarization tensor [34, 24]. To obtain the
isotropic conductivity of a moving plane we use the approach suggested in Ref.
[24] for graphene and take the limit, where the velocity Fermi and mass gap
tend to be zero. In this limit, the conductivity tensor becomes isotropic in
the co-moving frame of the plane. The conductivity in the laboratory is not
diagonal and depends on the velocity.

Throughout the paper the units $\hbar=c=1$ are used.

##  2 The Casimir energy of a moving plane

We use the approach for the Casimir effect of lateral moving graphene
developed in Ref. [24]. We consider two parallel conductive planes with
isotropic conductivities and inter-plane distance $a$. The first plane is at
rest in the laboratory frame and the second one is in a lateral motion with
velocity $\bm{v}$. The fluctuating electric field produces a current in the
plane according to the Ohm law. This current contributes to the boundary
condition and ultimately changes the energy spectrum. The second conductive
and moving plane is described by the same Ohm law in its co-moving frame. To
solve the scattering problem in the laboratory frame we need for conductivity
of the second plane in the laboratory frame. This question in relation to the
Ohm law was considered in Ref. [31, 32] by using the linear response tensor
method [33]. A similar approach was used in Ref. [24] where the polarization
tensor plays the role of the linear response tensor.

In the framework of the scattering matrix approach [27] the Casimir energy
$\mathcal{E}$ and pressure $\mathcal{P}$ for real frequencies may be
represented in the following form

$\mathcal{E}=-\frac{1}{2\mathrm{i}}\ESINT@iintop\ilimits@\frac{\mathrm{d}^{2}k}{(2\pi)^{3}}\left(I_{-}-I_{+}\right),\
\mathcal{P}=\ESINT@iintop\ilimits@\frac{\mathrm{d}^{2}k}{(2\pi)^{3}}\left(J_{-}+J_{+}\right),$
(1)

where

$\displaystyle I_{\pm}$
$\displaystyle=\int_{k}^{\infty}\mathrm{d}\omega\ln\det\left[1-e^{\pm
2\mathrm{i}ak_{3}}\bm{\mathcal{R}}(\pm k_{3})\right],\allowdisplaybreaks$
$\displaystyle J_{\pm}$ $\displaystyle=\int_{k}^{\infty}\mathrm{d}\omega
k_{3}\frac{e^{\pm 2\mathrm{i}ak_{3}}(\operatorname{tr}\bm{\mathcal{R}}(\pm
k_{3})-2e^{\pm 2\mathrm{i}ak_{3}}\det\bm{\mathcal{R}}(\pm
k_{3}))}{\det\left[\mathbf{I}-e^{\pm 2\mathrm{i}ak_{3}}\bm{\mathcal{R}}(\pm
k_{3})\right]},$ (2)

$\bm{\mathcal{R}}(\pm k_{3})=\bm{r}^{\prime}_{1}(\pm k_{3})\bm{r}_{2}(\pm
k_{3})$, and $k_{3}=\sqrt{\omega^{2}-\bm{k}^{2}}$. These formulas take into
account the propagating waves, only, because $\omega\geq k$. The subscript
$1(2)$ means that all reflection matrices related to the rest (moving) plane.
The scattering matrix for each part of system has the following form

$\bm{\mathcal{S}}=\begin{pmatrix}\bm{r}&\bm{t}^{\prime}\\\
\bm{t}&\bm{r}^{\prime}\end{pmatrix},$

with corresponding index and argument. The reflection matrices for conductive
plane were found in Ref. [27]

$\bm{r}_{i}=\bm{r}^{\prime}_{i}=-\frac{\omega^{2}\bm{\eta}_{i}-\bm{k}\otimes(\bm{k}\bm{\eta}_{i})+\mathbf{I}\omega
k_{3}\det\bm{\eta}_{i}}{\omega^{2}\operatorname{tr}\bm{\eta}_{i}-\bm{k}\bm{k}\bm{\eta}_{i}+\omega
k_{3}(1+\det\bm{\eta}_{i})},$ (3)

where $\bm{\eta}_{i}=2\pi\bm{\sigma}_{i}$, and $\bm{\sigma}_{i}$ is the
conductivity tensor of the plane $i=1,2$ ($i=1$ is at rest and $i=2$ is moving
in laboratory frame). In domain $\omega<k$ there are evanescent and waveguide
modes [35], but as shown in this paper, by rotation of the contour of
integration, the contribution of these modes are cancelled with energy of
boundary states of corresponding modes.

After the rotation of the contour to the imaginary axis, the two contributions
survived [24]. The first contribution with integration along the imaginary
frequency $\omega=\mathrm{i}\xi$ reads

$\mathcal{E}^{\perp}=\ESINT@iintop\ilimits@\frac{\mathrm{d}^{2}k}{2(2\pi)^{3}}\int_{-\infty}^{+\infty}\mathrm{d}\xi\ln\det\left[1-e^{-2ak_{E}}\bm{\mathcal{R}}(\mathrm{i}k_{E})\right],$
(4)

where $k_{E}=\sqrt{\xi^{2}+k^{2}}$. The corresponding force is perpendicular
to plane’s as usual Casimir force. This expression may be simplified by using
eigenvalues of matrix $\bm{\mathcal{R}}$ and represented as a sum of TE and TM
contributions. The matrices $\bm{r}^{\prime}_{1}$ and $\bm{r}_{2}$ are not
commutate and therefore the eigenvalues of $\bm{\mathcal{R}}$ are not a
product of the eigenvalues of $\bm{r}^{\prime}_{1}$ and $\bm{r}_{2}$. The
eigenvalues of $\bm{\mathcal{R}}$ may be found in closed but complicated forms
[24], which correspond to the contributions of TM and TE modes separately.
Instead of this approach, we use the expression for energy in the form [27]
directly via conductivity matrices $\bm{\eta}_{i}$

$\displaystyle\mathcal{E}^{\perp}$
$\displaystyle=\int\frac{\mathrm{d}^{2}k}{2(2\pi)^{3}}\int_{-\infty}^{\infty}\mathrm{d}\xi\ln\left(1+e^{-4ak_{E}}\frac{\xi^{2}k_{E}^{2}}{b_{1}b_{2}}\det\bm{\eta}_{1}\det\bm{\eta}_{2}\right.\allowdisplaybreaks$
$\displaystyle-\left.e^{-2ak_{E}}\left[\frac{\xi^{2}k_{E}^{2}}{b_{1}b_{2}}\left[(1-\det\bm{\eta}_{1})(1-\det\bm{\eta}_{2})+\det(\bm{\eta}_{1}-\bm{\eta}_{2})\right]-\frac{\xi
k_{E}}{b_{1}}-\frac{\xi k_{E}}{b_{2}}+1\right]\right),$ (5)

where
$b_{i}=\xi^{2}\operatorname{tr}\bm{\eta}_{i}+(\bm{k}\bm{k}\bm{\eta}_{i})+\xi
k_{E}\bigl{(}1+\det\bm{\eta}_{i}\bigr{)}$. The second contribution
$\mathcal{E}_{\parallel}$ gives a contribution to the force along planes, the
Casimir friction.

In Ref. [24], the normal force was considered for moving graphene with the
following conductivity in a co-moving frame

$\bm{\eta}_{1}=\eta_{\mathsf{gr}}\frac{\tilde{k}}{\omega}\left(\mathbf{I}+v_{F}^{2}\frac{\bm{k}\otimes\bm{k}}{\tilde{k}^{2}}\right)\widetilde{\Phi}\left(\frac{\tilde{k}}{2m}\right),$
(6)

where $\eta_{\mathsf{gr}}=2\pi\sigma_{\mathsf{gr}}=\pi e^{2}/2$, $v_{F}$ is
the Fermi velocity, and

$\widetilde{\Phi}(y)=\frac{2\mathrm{i}}{\pi
y}\left\\{1-\frac{y^{2}+1}{y}\operatorname{arctanh}y\right\\},\
\tilde{k}=\sqrt{\omega^{2}-v_{F}^{2}\bm{k}^{2}}.$ (7)

The general structure of the moving plane’s conductivity tensor in the
laboratory frame has the following form [24]

$\bm{\eta}_{2}=i_{1}\mathbf{I}+i_{2}\bm{k}\otimes\bm{k}+i_{3}(\bm{k}\otimes\bm{v}+\bm{v}\otimes\bm{k}),$
(8)

where

$\displaystyle i_{1}$
$\displaystyle=\frac{\eta_{\mathsf{gr}}\widetilde{\Phi}^{\prime}}{\omega\bm{k}^{2}\tilde{k}^{\prime}}\left(\bm{k}^{2}\tilde{k}^{\prime
2}+\frac{1-v_{F}^{2}}{1-v^{2}}((\bm{k}\bm{v})^{2}-\bm{k}^{2}\bm{v}^{2})k_{3}^{2}\right),\allowdisplaybreaks$
$\displaystyle i_{2}$
$\displaystyle=\frac{\eta_{\mathsf{gr}}\widetilde{\Phi}^{\prime}}{\omega\bm{k}^{2}\tilde{k}^{\prime}}\left(\bm{k}^{2}v_{F}^{2}+\frac{1-v_{F}^{2}}{1-v^{2}}\bm{v}^{2}k_{3}^{2}\right),\allowdisplaybreaks$
$\displaystyle i_{3}$
$\displaystyle=\frac{\eta_{\mathsf{gr}}\widetilde{\Phi}^{\prime}}{\bm{k}^{2}\tilde{k}^{\prime}}\frac{1-v_{F}^{2}}{1-v^{2}}\left(\bm{k}^{2}-\omega(\bm{k}\bm{v})\right),$
(9)

and

$\tilde{k}^{\prime}=\sqrt{\tilde{k}^{2}+\frac{1-v_{F}^{2}}{1-v^{2}}(\omega^{2}\bm{v}^{2}+(\bm{k}\bm{v})^{2}-2\omega(\bm{k}\bm{v}))},$
(10)

is the Lorentz transformation of $\tilde{k}$ (7).

##  3 The case of an isotropic conductivity

The isotropic conductivity case may be obtained by the formal limits $v_{F}\to
0$ and $\widetilde{\Phi}\to 1\ (m\to 0)$ in Eq. (6) and changing
$\eta_{\mathsf{gr}}$ to conductivity of corresponding plane $\eta_{i}$. After
these limits the conductivity tensor for the plane at rest becomes diagonal
$\bm{\eta}_{1}=\eta_{1}\mathbf{I}$. Taking the limits we obtain
$\tilde{k}^{\prime}=\gamma\omega_{v}$, where $\omega_{v}=\omega-\bm{k}\bm{v}$,
$\gamma=1/\sqrt{1-v^{2}}$ is relativistic factor, and

$\bm{\eta}_{2}=i^{\prime}_{1}\mathbf{I}+i^{\prime}_{2}\bm{k}\otimes\bm{k}+i^{\prime}_{3}(\bm{k}\otimes\bm{v}+\bm{v}\otimes\bm{k}),$
(11a) where
$i^{\prime}_{1}=\frac{\eta_{2}\gamma\left(\bm{k}^{2}\omega_{v}^{2}+((\bm{k}\bm{v})^{2}-\bm{k}^{2}\bm{v}^{2})k_{3}^{2}\right)}{\bm{k}^{2}\omega\omega_{v}},\
i^{\prime}_{2}=\frac{\eta_{2}\gamma\bm{v}^{2}k_{3}^{2}}{\bm{k}^{2}\omega\omega_{v}},\
i^{\prime}_{3}=\frac{\eta_{2}\gamma}{\bm{k}^{2}\omega_{v}}(\bm{k}^{2}-\omega(\bm{k}\bm{v})).$
(11b)

The $\omega_{v}\gamma$ is frequency of photon in laboratory which was emitted
in co-moving frame.

The straightforward calculations at the imaginary axis $\omega=\mathrm{i}\xi$
give

$\displaystyle\ b_{1}$
$\displaystyle=\left(\eta_{1}\xi+k_{E}\right)\left(\xi+\eta_{1}k_{E}\right),\
b_{2}=\frac{\xi}{\gamma\xi_{v}}\left(\eta_{2}\gamma\xi_{v}+k_{E}\right)\left(\gamma\xi_{v}+\eta_{2}k_{E}\right),\allowdisplaybreaks$
$\displaystyle\det(\bm{\eta}_{1}-\bm{\eta}_{2})$
$\displaystyle=(\eta_{1}-\eta_{2})^{2}+\eta_{1}\eta_{2}\left(\gamma
v^{2}\frac{k_{E}^{2}}{\xi\xi_{v}}+2(1-\gamma)\right),\
\det\bm{\eta}_{i}=\eta_{i}^{2},$ (12)

where $\xi_{v}=\xi+\mathrm{i}\bm{k}\bm{v}$. Then we use the polar coordinates
for $\bm{k}$ in Eq. (5), $\bm{k}\bm{v}=kv\cos\varphi$ and transform
coordinates of plane $k\in[0,\infty)$, $\xi\in(-\infty,\infty)$ to polar
coordinates $\xi=k_{E}\cos\theta,k=k_{E}\sin\theta$ where $\theta\in[0,\pi]$.
After these changes, the dependence of $k_{E}$ is survived in exponents only.
By changing the variable $ak_{E}=y$ we observe that the energy depends on the
inter-plane distance as $1/a^{3}$ for constant conductivities, as expected
[36]. Thus, the energy has the following form ($x=\cos\theta$)

$\mathcal{E}^{\perp}_{1,2}=\mathrm{Re}\int_{0}^{\infty}\frac{y^{2}\mathrm{d}y}{(2\pi
a)^{3}}\int_{0}^{1}\mathrm{d}x\int_{0}^{\pi}\mathrm{d}\varphi E_{1,2},\
\mathcal{P}^{\perp}_{1,2}=\frac{3}{a}\mathcal{E}^{\perp}_{1,2},$ (13)

where

$\displaystyle E_{1,2}$
$\displaystyle=\ln\left(1+e^{-4y}\frac{x^{2}\eta_{1}^{2}\eta_{2}^{2}}{\beta_{1}\beta_{2}}\right.\allowdisplaybreaks$
$\displaystyle\left.-e^{-2y}\left(\frac{x^{2}}{\beta_{1}\beta_{2}}\left[(1-\eta_{1}\eta_{2})^{2}+\eta_{1}\eta_{2}\left(\gamma
v^{2}\frac{k_{E}^{2}}{\xi\xi_{v}}+2(1-\gamma)\right)\right]-\frac{x}{\beta_{1}}-\frac{x}{\beta_{2}}+1\right)\right),$
(14)

where

$\beta_{1}=\left(\eta_{1}x+1\right)\left(x+\eta_{1}\right),\
\beta_{2}=\frac{x}{\gamma x_{v}}\left(\eta_{2}\gamma
x_{v}+1\right)\left(\gamma x_{v}+\eta_{2}\right),\
x_{v}=x-\mathrm{i}v\sqrt{1-x^{2}}\cos\varphi.$ (15)

If the first plane (at rest) is a perfect conductor we tend
$\eta_{1}\to\infty$ and obtain

$E_{\mathrm{id},2}=\ln\left(1+e^{-4y}\frac{x\eta_{2}^{2}}{\beta_{2}}-e^{-2y}\left(\frac{x\eta_{2}^{2}}{\beta_{2}}-\frac{x}{\beta_{2}}+1\right)\right),$
(16)

and for two ideal planes

$E_{\mathrm{id},\mathrm{id}}=2\ln\left(1+e^{-2y}\right),$ (17)

the energy does not depend on the velocity.

Without a relative movement, $\bm{v}=0$, we return to the results obtained in
Ref. [36]:

$E=\ln\left(1-e^{-2y}\frac{\eta_{1}\eta_{2}}{(\eta_{1}+x)(\eta_{2}+x)}\right)+\ln\left(1-e^{-2y}\frac{\eta_{1}\eta_{2}x^{2}}{(x\eta_{1}+1)(x\eta_{2}+1)}\right)=E_{\mathsf{tm}}+E_{\mathsf{te}},$
(18)

the sum of TM and TE contributions.

Figure 1: The plots of the velocity correction
$\frac{\Delta_{v}\mathcal{E}_{\perp}}{\mathcal{E}_{0}}=\frac{\mathcal{E}_{\perp}-\mathcal{E}_{0}}{\mathcal{E}_{0}}$
for $\eta=0.01,0.1$. For small velocity and conductivity $v\ll\eta\ll 1$ (left
and middle panels) we have quadratic correction according with (19). The solid
lines are for two planes with the same conductivity (14) and dashed lines are
for system with a plane at rest with a perfect conductivity (16).

Let us consider the constant conductivities case $\eta_{1}=\eta_{2}=\eta$. For
small velocity and conductivity, and $v\ll\eta\ll 1$ we obtain from Eq. (13):

$\mathcal{E}^{\perp}_{i,k}\approx
e_{i,k}\frac{\eta}{a^{3}}\left(\frac{v}{\eta}\right)^{2},\mathcal{P}^{\perp}_{i,k}\approx
p_{i,k}\frac{\eta}{a^{4}}\left(\frac{v}{\eta}\right)^{2},$ (19)

where

$e_{1,2}=4.8\cdot 10^{-4},\ e_{\mathrm{id},2}=6.1\cdot 10^{-4},\
p_{1,2}=1.4\cdot 10^{-3},\ p_{\mathrm{id},2}=1.8\cdot 10^{-3}.$ (20)

Numerical evaluations of (13) are shown in Fig. 1 for two systems: $(1,2)$ –
two conductive planes with constant conductivity $\eta$ (solid lines) and
$(\mathrm{id},2)$ – the first plane at rest is a perfect metal (dashed lines).
We calculated velocity correction to the energy:
$\Delta_{v}\mathcal{E}_{\perp}/\mathcal{E}_{0}=(\mathcal{E}_{\perp}-\mathcal{E}_{0})/\mathcal{E}_{0}$,
where $\mathcal{E}_{0}$ is energy without relative movement. The velocity
correction is negative for both systems and it reveals quadratic behaviour
$\sim v^{2}$ for small velocity $v\ll\eta\ll 1$. The absolute value of
correction is greater for the first system.

The above derivation is applicable for isotropic but frequency-dependent
conductivity, $\eta=\eta(\omega)$. Let us consider the simple case of Drude
conductivity with

$\eta_{1}=\frac{\eta\Gamma}{\Gamma+\xi}\ ,\
\eta_{2}=\frac{\eta\Gamma}{\Gamma+\gamma\xi_{v}},$ (21)

and parameters are taken for graphene, $\Gamma=6.365$ eV and
$\eta=\eta_{\mathsf{gr}}=e^{2}/4$ [37]. After changing the integrand variables
as were made above, the Casimir energy acquires an additional dependence on
the inter-plane distance through conductivity:

$\eta_{1}=\frac{\eta_{\mathsf{gr}}(a\Gamma)}{(a\Gamma)+yx}\ ,\
\eta_{2}=\frac{\eta_{\mathsf{gr}}(a\Gamma)}{(a\Gamma)+\gamma yx_{v}}.$ (22)

For $a=100$ nm and $\Gamma=6.365$ eV, one has, $a\Gamma=3.225$. For
$a\Gamma\gg 1$ the conductivities $\eta_{1}=\eta_{2}=\eta_{\mathsf{gr}}$, as
should be the case – the constant conductivity model is valid for large inter-
plane distances.

Figure 2: The plots of velocity correction
$\frac{\Delta_{v}\mathcal{E}_{\perp}}{\mathcal{E}_{0}}=\frac{\mathcal{E}_{\perp}-\mathcal{E}_{0}}{\mathcal{E}_{0}}$
for Drude model ($a=10,100$ nm) and constant conductivity with
$\eta=\eta_{\mathsf{gr}}$. Left panel: the region of small velocity $v\ll 1$.
Right panel: the whole range of $v$. The solid lines are for two planes with
the same conductivity and dashed lines are for system where a plane at rest is
a perfect conductor.

The numerical evaluations of Casimir energy are plotted in Fig. 2. We observe
that the greater the inter-plane distance the closer energy is to the case of
constant conductivity (blue lines) as should be the case.

Let us compare the obtained results with the case of graphene analysed in Ref.
[24]. For two graphene sheets the energy as a function of velocity at the
beginning is positive and then becomes negative with maximum for
$v_{c}=v_{F}+(ma)/2$. The region with positive energy disappears for
$m=v_{F}=0$. The case of two planes with constant conductivity reveals a
negative value of energy for all values of velocity. The velocity correction
$\sim v^{2}$ for both cases. Concerning inter-plane distance we have different
pictures. For two planes with constant conductivity, the energy is $\sim
1/a^{3}$ for any distance while for graphenes it has this dependence for large
distances only. It is expected because the constant conductivity model for
graphene is valid for large inter-plane distances. For the case of a system
perfect conductor/graphene the energy is zero up to a specific velocity, while
for the above-considered case, we have quadratic behaviour at the beginning. A
similar conclusion is valid for the case of the Drude model of conductivity
(21). For large inter-plane distances, both models are very close. For small
distances, the weak dependents on the distance takes place.

##  4 Conclusion

In this paper we considered the normal (perpendicular to the planes) Casimir
force for two conductive planes with an isotropic conductivity that is
laterally moving with relative velocity $v$. The main problem is connected
with calculating the conductivity of a moving plane in a laboratory frame. In
a co-moving frame, the isotropic conductivity is a coefficient in the Ohm law,
$\bm{J}^{\prime}=\sigma^{\prime}\bm{E}^{\prime}$, where $\bm{E}^{\prime}$ and
$\bm{J}^{\prime}$ are fluctuations of the electric field and corresponding
density of current. Transformation of this relation to the laboratory frame
(where the first plane is at the rest) is not trivial. The simple way [31, 32]
to solve this problem is to start from the linear relation between the density
of current and electromagnetic vector potential,
$J^{\mu}=\Pi^{\mu}_{\nu}A^{\nu}$ which is usual for plasma physics [33]. The
$\Pi^{\mu}_{\nu}$ is the linear response tensor. The same approach was used in
Ref. [34] where the polarization tensor played the role of the linear response
tensor $\Pi^{\mu}_{\nu}$.

The transformation of the conductivity has no simple form [32] even for a
constant conductivity case. A similar approach was used in Ref. [24] where the
linear relation for current and electromagnetic potential, and invariance of
the boundary condition were used. To obtain the case of isotropic conductivity
we used expressions obtained for graphene in which we take limits $v_{F}\to 0$
for Fermi velocity and $m\to 0$ for mass gap. With these limits the
conductivity tensor in the co-moving frame becomes diagonal. In the laboratory
frame, it has the form (11). Having these tensors, the Casimir energy may be
calculated by expression (5), firstly obtained in Ref. [27].

Analysis of expressions obtained for two conductive planes (14), and for
system (perfect conductivity)/(constant conductivity) is difficult due to two
small dimensionless parameters: velocity of plane $v$ and conductivity
$\eta=2\pi\sigma$ (dimensionless for 2D systems). For the case $v\ll\eta\ll 1$
the energy $\sim\eta(v/\eta)^{2}$ (19). The quadratic dependence is usual for
normal force and different directions of motion [7, 8, 10, 24]. The inter-
plane distance dependence of energy is $1/a^{3}$ for any distance as usual for
constant conductivity case [36] because the constant conductivity model is
valid for large distances where the Casimir regime is satisfied. The Drude
model of conductivity reveals close behaviour of the system with a weak
dependence on the inter-plane distances (see Fig. 2).

## Acknowledgments

NK was supported in part by the grants 2022/08771-5, 2021/10128-0 of São Paulo
Research Foundation (FAPESP).

## References

  * Casimir [1948] H. B. G. Casimir, Kon. Ned. Akad. Wetensch. Proc. 51, 793 (1948).
  * Bordag  _et al._ [2009a] M. Bordag, G. L. Klimchitskaya, U. Mohideen, and V. M. Mostepanenko,  _Advances in the Casimir Effect_ (Oxford University Press, 2009) p. 768.
  * Woods  _et al._ [2016] L. M. Woods, D. A. R. Dalvit, A. Tkatchenko, P. Rodriguez-Lopez, A. W. Rodriguez, and R. Podgornik, Rev. Mod. Phys. 88, 45003 (2016).
  * Dodonov [2020] V. Dodonov, Physics 2, 67 (2020).
  * Bordag  _et al._ [2001] M. Bordag, U. Mohideen, and V. Mostepanenko, Phys. Rep. 353, 1 (2001).
  * Khusnutdinov and Woods [2019] N. Khusnutdinov and L. M. Woods, JETP Lett. 110, 183 (2019).
  * Bordag  _et al._ [1984] M. Bordag, G. Petrov, and D. Robaschik, Sov. J. Nucl. Phys. 39, 828 (1984).
  * Bordag  _et al._ [1986] M. Bordag, F. Dittes, and D. Robaschik, Sov. J. Nucl. Phys. 43, 1034 (1986).
  * Fulling and Davies [1976] S. Fulling and P. Davies, Proc. R. Soc. Lond. A 348, 393 (1976).
  * Maslovski [2011] S. I. Maslovski, Phys. Rev. A 84, 22506 (2011).
  * Polevoi [1990] V. Polevoi, Sov. Phys. JETP 71, 1119 (1990).
  * Dedkov and Kyasov [2018] G. Dedkov and A. Kyasov, Chinese Journal of Physics 56, 3002 (2018).
  * Mkrtchian [1995] V. Mkrtchian, Phys. Lett. A 207, 299 (1995).
  * Pendry [1997] J. B. Pendry, J. Phys.: Condens. Matter 9, 10301 (1997).
  * Volokitin and Persson [2007] A. I. Volokitin and B. N. J. Persson, Rev. Mod. Phys. 79, 1291 (2007).
  * Philbin and Leonhardt [2009] T. G. Philbin and U. Leonhardt, New J. Phys. 11, 033035 (2009).
  * Pendry [2010] J. B. Pendry, New J. Phys. 12, 033028 (2010).
  * Maghrebi  _et al._ [2013a] M. F. Maghrebi, R. Golestanian, and M. Kardar, Phys. Rev. D 87, 025016 (2013a).
  * Maghrebi  _et al._ [2013b] M. F. Maghrebi, R. Golestanian, and M. Kardar, Phys. Rev. A 88, 042509 (2013b).
  * Farias  _et al._ [2017] M. B. Farias, C. D. Fosco, F. C. Lombardo, and F. D. Mazzitelli, Phys. Rev. D 95, 065012 (2017).
  * Farias  _et al._ [2018] M. B. Farias, W. J. M. Kort-Kamp, and D. A. R. Dalvit, Phys. Rev. B 97, 161407 (2018).
  * Brevik  _et al._ [2022] I. Brevik, B. Shapiro, and M. G. Silveirinha, Int. J. Mod. Phys. A 37, 22410123 (2022).
  * Guo  _et al._ [2023] X. Guo, K. A. Milton, G. Kennedy, and N. Pourtolami, Phys. Rev. A 107, 062812 (2023).
  * Antezza  _et al._ [2023] M. Antezza, N. Emelianova, and N. Khusnutdinov, Casimir-Lifshitz force for moving graphene (2023).
  * Levin and Rytov [1967] M. Levin and S. Rytov,  _Theory of Equilibrium Thermal Fluctuations in Electrodynamics_ (Nauka, Moscow, in Russian, 1967).
  * Zel’dovich [1971] Ya.B. Zel’dovich, JETP Lett. 14, 270 (1971).
  * Fialkovsky  _et al._ [2018] I. Fialkovsky, N. Khusnutdinov, and D. Vassilevich, Phys. Rev. B 97, 165432 (2018).
  * Van Kampen  _et al._ [1968] N. Van Kampen, B. Nijboer, and K. Schram, Physics Letters A 26, 307 (1968).
  * Gerlach [1971] E. Gerlach, Phys. Rev. B 4, 393 (1971).
  * Silveirinha [2014] M. G. Silveirinha, Phys. Rev. X 4, 031013 (2014).
  * Starke and Schober [2015] R. Starke and G. Schober, Photonics Nanostructures: Fundam. Appl. 14, 1 (2015).
  * Starke and Schober [2016] R. Starke and G. A. H. Schober, Int. J. Mod. Phys. D 25, 1640010 (2016).
  * Melrose [2008] D. B. Melrose,  _Quantum Plasmadynamics: Unmagnetized Plasmas_ , Lecture Notes in Physics No. 735 (Springer, New York, 2008).
  * Bordag  _et al._ [2009b] M. Bordag, I. V. Fialkovsky, D. M. Gitman, and D. V. Vassilevich, Phys. Rev. B 80, 245406 (2009b).
  * Bordag [2012] M. Bordag, Phys. Rev. D 85, 025005 (2012).
  * Khusnutdinov  _et al._ [2014] N. Khusnutdinov, D. Drosdoff, and L. M. Woods, Phys. Rev. D 89, 085033 (2014).
  * Khusnutdinov  _et al._ [2015] N. Khusnutdinov, R. Kashapov, and L. M. Woods, Phys. Rev. D 92, 045002 (2015).