PDE-constrained optimal control of a leader-follower opinion formation model

TL;DR

This paper develops a PDE-constrained optimal control framework for leader-follower opinion formation models, deriving associated equations and solving them numerically to analyze different interaction scenarios.

Contribution

It introduces a novel PDE-constrained optimal control approach for leader-follower opinion dynamics, including derivation of optimality conditions and a numerical solution method.

Findings

Effective numerical algorithms for opinion control.

Simulation results for various interaction types.

Insights into influence of control strategies.

Abstract

We consider the PDE-constrained optimal control of a leader-follower kinetic opinion formation model, with a Fokker-Planck-type system of partial differential equations as a state constraint. We derive the Boltzmann-type and Fokker-Planck-type systems of equations associated with the controlled leader-follower opinion formation model. In a function space setting we derive first-order optimality conditions associated with the PDE-constrained optimal control problem, yielding an optimality system of coupled nonlinear partial differential equations. We employ a gradient-type sweeping algorithm to numerically attack the optimality system obtained from the first-order optimality conditions. We present the results from a finite elements based simulation for different types of interactions and cost functionals.

Tables

Table 1. Table 1 : Default parameters used for the Fokker-Planck simulations.

Parameter	${\tau }_{LL}$	${\tau }_{FL}$	${\tau }_{FF}$	${\lambda }_L$	${\lambda }_F$	$r$	$\beta$
Value	0.2	2	0.2	0.05	0.05	0.5	0.05

Equations

$$\begin{array}{rl}{w}_{Li}^{\ast }& =w_{Li}+{\gamma }_L{P}_L(w_{Li},w_{Lj})(w_{Lj}-w_{Li})+{\eta }_{L}D(w_{Li})\\ & +\frac{{\gamma }_L}{2}u(w_{Li},w_{Lj},t)\\ w_{Lj}^{\ast }& =w_{Lj}+{\gamma }_L{P}_L(w_{Lj},w_{Li})(w_{Li}-w_{Lj})+{\tilde{\eta }}_{L}D(w_{Lj})\\ & +\frac{{\gamma }_L}{2}u(w_{Lj},w_{Li},t)\end{array}$$

$$\begin{array}{rl}{w}_{Li}^{\ast }& =w_{Li}\\ w_{Fj}^{\ast }& =w_{Fj}+{\gamma }_F\tilde{P}(w_{Fj},w_{Li})(w_{Li}-w_{Fj})+\tilde{\eta }D(w_{Fj})\end{array}$$

$$\begin{array}{rl}{w}_{Fi}^{\ast }& =w_{Fi}+{\gamma }_F{P}_F(w_{Fi},w_{Fj})(w_{Fj}-w_{Fi})+{\eta }_{F}D(w_{Fi})\\ w_{Fj}^{\ast }& =w_{Fj}+{\gamma }_F{P}_F(w_{Fj},w_{Fi})(w_{Fi}-w_{Fj})+{\tilde{\eta }}_{F}D(w_{Fj})\end{array}$$

$$D(x)=(1-x^2{)}^{\alpha },$$

$$M_{2+\delta }$$

$$\displaystyle\qquad\left\langle|\eta|^{\alpha}\right\rangle=\int_{I}|\eta|^{\alpha}\ \text{d}\mathbb{P}(\eta)<\infty,\forall\alpha\leq 2+\delta\bigg{\}},$$

$$\begin{split}\frac{\text{d}}{\text{d}t}\left(\int_{\mathcal{I}}\phi(w)f_{L}(w,t)\ \text{d}w\right)=&\frac{1}{\tau_{LL}}\Big{\langle}\int_{\mathcal{I}^{2}}\big{[}\phi(w^{*})+\phi(v^{*})-\phi(w)-\phi(v)\big{]}\\ &\qquad\qquad\times f_{L}(w,t)f_{L}(v,t)\ \text{d}w\text{d}v\Big{\rangle}.\end{split}$$

$$\begin{split}\frac{\text{d}}{\text{d}t}\left(\int_{\mathcal{I}}\phi(w)f_{F}(w,t)\ \text{d}w\right)=&\frac{1}{\tau_{FF}}\Big{\langle}\int_{\mathcal{I}^{2}}\big{[}\phi(w^{*})+\phi(v^{*})-\phi(w)-\phi(v)\big{]}\\ &\qquad\qquad\times f_{F}(w,t)f_{F}(v,t)\ \text{d}w\text{d}v\Big{\rangle}\\ &+\frac{1}{\tau_{FL}}\Big{\langle}\int_{\mathcal{I}^{2}}\big{[}\phi(w^{*})-\phi(w)\big{]}\\ &\qquad\qquad\times f_{F}(w,t)f_{L}(v,t)\ \text{d}w\text{d}v\Big{\rangle}.\end{split}$$

$$\begin{array}{rl}\frac{\partial g_L}{\partial s}(w,s)=& \frac{\partial }{\partial w}\left(\left(\frac{1}{{\tau }_{LL}}\mathcal{K}[g_L](w,s)+\frac{1}{2{\tau }_{LL}}u(w,s)\right)g_L(w,s)\right)\\ & +\frac{{\lambda }_L}{2{\tau }_{LL}}\frac{{\partial }^2}{\partial w^2}(D^2(w)g_L(w,s)),\end{array}$$

$$\mathcal{K}[g_L](w,s):={\int }_{\mathcal{I}}{P}_L(w,v)(w-v)g_L(v,s)\text{d}v,$$

$$\begin{array}{rl}\frac{\partial g_F}{\partial s}(w,s)=& {\alpha }_{LF}\frac{\partial }{\partial w}\left(\left(\frac{1}{2{\tau }_{FL}}\mathcal{M}[g_L](w,s)+\frac{1}{{\tau }_{FF}}\mathcal{N}[g_F](w,s)\right)g_F(w,s)\right)\\ & +\left(\frac{{\lambda }_F}{4{\tau }_{LF}}+\frac{{\lambda }_F}{2{\tau }_{FF}}\right)\frac{{\partial }^2}{\partial w^2}(D^2(w)g_F(w,s)),\end{array}$$

$$\mathcal{M}[g_L](w,s):=$$

$$\mathcal{N}[g_F](w,s):=$$

$$\frac{\partial\mathbf{g}}{\partial t}=\frac{\partial}{\partial w}\bigg{(}(A[\mathbf{g}]+B[u])\mathbf{g}+\frac{\partial}{\partial w}(C\mathbf{g})\bigg{)},$$

$$\mathbf{g}$$

$$A[\mathbf{g}]$$

$$B[u]$$

$$C$$

$$\mathbf{g}(0,w)$$

$$\begin{array}{r}\underset{u\in U}{\min }J(\mathbf{g},u),\end{array}$$

$$\frac{\partial\mathbf{g}}{\partial t}=\frac{\partial}{\partial w}\bigg{(}(A[\mathbf{g}]+B[u])\mathbf{g}+\frac{\partial}{\partial w}(C\mathbf{g})\bigg{)},$$

$$\mathbf{g}(0,w)={\mathbf{g}}^0(w),$$

$$J(\mathbf{g},u)=\frac{1}{2}{\int }_0^T{\int }_{\mathcal{I}}\left(|w-w_{dF}{|}^2+\beta |u{|}^2\right)g_F\text{d}w\text{d}t,$$

$$U=\big{\{}||u||_{L^{2}(0,T;L^{\infty}(\mathcal{I}))}\leq M:u\in L^{2}(0,T;L^{\infty}(\mathcal{I}))\big{\}},$$

$$\displaystyle\int_{0}^{T}\Big{\langle}\frac{\partial g_{L}}{\partial t},\varphi\Big{\rangle}_{H^{-1},H^{1}}dt-$$

$$+\frac{{\lambda }_L}{2{\tau }_{LL}}\frac{\partial }{\partial w}(D^2(w)g_L)\frac{\partial \phi }{\partial w}dwdt=0,$$

$$\displaystyle\int_{0}^{T}\Big{\langle}\frac{\partial g_{F}}{\partial t},\varphi\Big{\rangle}_{H^{-1},H^{1}}dt-$$

$$+\left(\frac{{\lambda }_F}{4{\tau }_{LF}}+\frac{{\lambda }_F}{2{\tau }_{FF}}\right)\frac{\partial }{\partial w}(D^2(w)g_F(w,s))\frac{\partial \phi }{\partial w}dwdt=0.$$

$$\underset{j\to \infty }{\lim }J({\mathbf{g}}^j,u^j)=\underset{u\in U}{\inf }J(\mathbf{g},u).$$

$$g_L^{j_k},g_F^{j_k}\to {\bar{g}}_L,{\bar{g}}_F\text{in}{L}^2(\mathcal{I}\times (0,T))\text{and}{u}^{j_k}\stackrel{\ast }{\rightharpoonup }\bar{u}\text{in}{L}^2(0,T;L^{\infty }(\mathcal{I}))\text{as}k\to \infty .$$

$$\underset{u\in U}{\inf }J(\mathbf{g},u)=\underset{j\to \infty }{\lim }J({\mathbf{g}}^j,u^j)=\underset{k\to \infty }{\mathrm{liminf}}J({\mathbf{g}}^{j_k},u^{j_k})\ge J(\bar{\mathbf{g}},\bar{u}),$$

$$\mathcal{L}(\mathbf{g},u,\mathbf{p})=J(\mathbf{g},u)-\left(\frac{\partial \mathbf{g}}{\partial s},\mathbf{p}\right)+\left(\frac{\partial }{\partial w}((A[\mathbf{g}]+B[u])\mathbf{g}),\mathbf{p}\right)-\left(\frac{\partial }{\partial w}(C\mathbf{g}),\frac{\partial \mathbf{p}}{\partial w}\right),$$

$${\nabla }_{\mathbf{p}}\mathcal{L}(\bar{\mathbf{g}},\bar{u},\bar{\mathbf{p}})\varphi$$

$${\nabla }_{\mathbf{g}}\mathcal{L}(\bar{\mathbf{g}},\bar{u},\bar{\mathbf{p}})\psi$$

$${\nabla }_u\mathcal{L}(\bar{\mathbf{g}},\bar{u},\bar{\mathbf{p}})v$$

$$\begin{array}{r}{\nabla }_{\mathbf{p}}\mathcal{L}\varphi =\left(\frac{\partial \mathbf{g}}{\partial s},\varphi \right)-\left(\frac{\partial }{\partial w}((A[\mathbf{g}]+B[u])\mathbf{g}),\varphi \right)+\left(\frac{\partial }{\partial w}(C\mathbf{g}),\frac{\partial \varphi }{\partial w}\right)=0,\end{array}$$

Full text

MSC Classification]49M41, 91D30, 35Q93

1 Introduction

We study the PDE-constrained optimal control of a leader-follower kinetic opinion formation model from [1], which features two species, a species of strong opinion leaders and a species of followers, building on the works of Fornasier and Solombrino [2] and Albi et al. [3]. In [3] the authors consider the mean field optimal control of a partial differential equation of continuity type, as they arise in opinion formation problems for a large population in the mean-field limit for a single species. They derive first-order optimality conditions and propose a sub-optimal, but efficient numerical optimal control strategy.

The introduction of opinion formation leader-follower models in the literature in the past two decades [4, 1, 5] has naturally sparked curiosity on how leaders can optimise their behaviour, e.g. position themselves optimally to attract a mass of followers. Another question that can be asked is how leaders can optimally behave to steer a follower density to a desired target distribution. Examples in the real world spring to mind, with political parties, their leaders and their spin doctors, trying to attract voters, e.g. by pushing certain topics where differentiation from other parties is easier, or to establish them as more visibly different alternative to other parties.

PDE constrained optimal control in an infinite dimensional setting has been dynamically developing as a research area since the late 1990s and 2000s years (confer [6, 7] for an overview and references), moving from prototype nonlinear PDE to systems of nonlinear PDE, e.g. from fluid dynamics, and continues to thrive. For the nonlinear, nonlocal partial differential equations arising in opinion formation, many works have favoured numerically less expensive approaches. One reason is that the nonlocal terms in these equations make their numerical solution and even more so their numerical optimisation very challenging. A popular class of algorithms is known as model predictive control, which is a powerful, yet sub-optimal approach to attack the complex control problem [8, 9] by solving iteratively small finite-horizon optimal control problems. Similar many-particle models and associated control approaches are present in flocking [10, 11], crowd dynamics [12, 13, 14], epidemic dynamics [15, 16, 17], opinion formation on graphs [18], traffic flow [19] and global optimisation [20, 21].

In this paper, starting from controlled interactions in Section 2, we derive the Boltzmann-type and Fokker-Planck-type systems of partial differential equations associated with the controlled leader-follower opinion formation model in Section 3. Section 4 is devoted to PDE-constrained optimal control. In an infinite-dimensional function space setting we derive first-order optimality conditions associated with the PDE-constrained optimal control problem with the Fokker-Planck-type system as a state constraint, yielding an optimality system of coupled nonlinear partial differential equations. In Section 5 we discuss the numerical algorithms. We employ a gradient-type sweeping algorithm [22] to numerically attack the optimality system obtained from the first-order optimality conditions. We present the results from a finite elements based simulation for different types of interactions and cost functionals in Section 6.

2 Multi-agent opinion formation model

We consider a large number of interacting individuals, with each individual’s opinion in an interval $\mathcal{I}:=[-1,1]$. The overall population consists of two species, an opinion leader species and a follower species. Similar as in [1], members of either species can interact with each other through binary interactions, with three possible cases: (a) the interactions where a member of the leader species interacts with another member of the leader species (an L and L interaction); (b) a member of the leader species interacts with a member of the follower species (an L and F interaction); and (c) a member of a follower species interacts with another member of a follower species (an F and F interaction). We assume that the leader species is formed by strong opinion leaders, as in [1], that is assertive individuals being able to influence followers but withstanding being influenced themselves. To that end, we have that an individual of the leader species will have their opinion unchanged by an interaction with an individual from a follower species, but interact normally with another leader. Unlike [1] we allow a control function $u$ to modify the leader interactions. The interactions are then as follows,

• •

L-L interactions:

[TABLE]

• •

L-F interactions:

[TABLE]

• •

F-F interactions:

[TABLE]

where $w_{Li}^{*}$, $w_{Fi}^{*}$ are the post-interaction opinions of the individuals.

Note here that we add the control term only to the leader interaction (1) which represents the effect of a particular strategy that involves changing the leader’s apparent opinion in order to minimise a given cost functional $J$. This cost functional may be as simple as trying to maximise a follower base – as is likely in the case of some form of marketing strategy – or, in a paradigm where multiple leader species exist, to minimise the follower base of the other leader species. Since we are only controlling the leader interactions we expect the action of the control function $u$ on the follower species to be only through the interactions with leader species.

In the terms with the form $\gamma_{\nu}P_{\nu}(w,v)(v-w)$ and $\gamma_{\nu}\tilde{P}(w,v)(v-w)$, we are representing a compromise between the opinion of the first and second individuals. In the prevalent literature, the $P$ functions in these interactions has been taken as either a bounded confidence model – $P$ functions are taken as localisation functions such as a characteristic function or a hyperbolic tangent function [23, 5, 1] – using the assumption that individuals will only interact with other individuals with similar opinions, or the Sznajd model [24]– $P$ functions penalise extremal opinions – where it is assumed that more extreme opinions interact little with any other opinions. A more in-depth discussion of these choices for compromise functions is given below in Section 5. More realistic choices for this compromise function while also taking into account the different “demographics" of the population are proposed in [25, 26].

The constants $\gamma_{L}$ and $\gamma_{F}$ in these terms govern the rate of compromise over time. We note here that we choose $P$ such that the range of $P$ is in $[0,1]$ and we choose $\gamma_{L}$ and $\gamma_{F}$ to be in $[0,\tfrac{1}{2})$ to ensure post-interaction opinions remain in $\mathcal{I}$ (see [1]).

The terms of the form $\eta_{\nu}D(w)$, represent the influence of other sources, which in our case act randomly. This includes things like, media controlled by outside sources, the individual thinking about their opinion. Here we are making the assumption that these outside influences are unbiased. Here $\eta_{L}$ and $\eta_{F}$ represent the random elements. The function $D$ is chosen in such a way as to ensure that the diffusion term is sufficiently small near the boundaries – and should not be confused with a differential operator. This is to ensure that the post-interaction opinions of the individuals $w^{*}$, $v^{*}$ remain within $\mathcal{I}$. A usual choice for this function in the literature [23, 5, 1] is,

[TABLE]

for some $\alpha>0$. We use this choice of $D$ with $\alpha=2$ throughout our numerical experiments.

The random variables $\eta_{L}$ and $\eta_{F}$ are chosen such that their probability density function $\Theta$ from the set of probability measures,

[TABLE]

for some $\delta>0$. We can construct this density a random variable $Y$ with the properties, $\left\langle Y\right\rangle=\int_{I}Y\ \text{d}\Theta(Y)=0$, $\left\langle Y^{2}\right\rangle=\int_{I}Y^{2}\ \text{d}\Theta(Y)=1$, then we choose $\eta_{\nu}$ with the property that $\left\langle|\eta_{\nu}|^{p}\right\rangle=\left\langle|Y\sigma_{\nu}|^{p}\right\rangle$ for $\nu=L,F$, and $p\in[0,2+\delta]$, with $0<\sigma_{\nu}<1$. Here $\sigma_{\nu}^{2}$ is the variance of the random variables $\eta_{\nu}$.

3 Partial differential equations

In this section, we discuss the system of Boltzmann-type equations that results from the interactions (1)–(3), and then discuss the system of Fokker-Planck-type equations, which are obtained from the quasi-invariant limit of the Boltzmann-type system. Throughout the rest of this paper, we make the simplifying assumption $u(w,v,t)=u(w,t)$.

3.1 Boltzmann-type system

Using standard methods of kinetic theory we can now derive the Boltzmann-type system associated with the binary interactions, for completeness this process is outlined in Appendix A. Let $\phi$ be test functions from $C_{c}^{2,\delta}(\mathcal{I};\real)$, that are continuous and compactly supported functions on $\mathcal{I}$.

The Boltzmann-type equation for the leader species $f_{L}(w,t)$ in weak form is given by

[TABLE]

For the follower species $f_{F}(w,t)$ we have the Boltzmann-type equation in weak form,

[TABLE]

In the above we denote by $\tau_{LL}$, $\tau_{FL}$ and $\tau_{FF}$ the relaxation times associated with the respective interactions.

3.2 Fokker-Planck-type system

In the so-called quasi-invariant limit a coupled system of non-linear, non-local Fokker-Planck equations is obtained from (4), (5) for this model. To this end time is rescaled setting $s=\gamma_{L}t$ and the transformed density functions $g_{L}(w,s)=f_{L}(w,t)$ are considered, taking the limit $\sigma,\gamma_{L}\rightarrow 0$ while keeping the value $\sigma^{2}/\gamma_{L}$ constant, see Appendix B for further details. We then obtain the Fokker-Planck equation for the leader species,

[TABLE]

where

[TABLE]

subject to no-flow boundary conditions.

Taking the same limit and setting $\alpha_{LF}=\gamma_{L}/\gamma_{F}$ (see Appendix B for further details), the follower Fokker-Planck-type equation is given by,

[TABLE]

where

[TABLE]

subject to no-flow boundary conditions.

The system (6) and (7) can be expressed as a matrix-vector equation, namely,

[TABLE]

with,

[TABLE]

4 PDE-constrained optimal control

We introduce the optimal control problem with the system of Fokker-Planck-type equations (8) as a PDE constraint. We derive the first-order optimality conditions using the formal Lagrangian Method.

4.1 Optimal control problem

We now define the control problem for this model. The optimal control of the model can come in three shapes: 1. microscopic – which optimises a cost functional with the microscopic interactions as the state; 2. mesoscopic – where the cost functional is minimised subject to the system of Boltzmann-type equations (4), (5); 3. macroscopic – which minimises the cost functional subject to the Fokker-Planck-type system. The hierarchy of these optimal control problems for a single species are discussed in [3]. Here, we are considering the optimal control of the system of Fokker-Planck-type equations arising for the two species leader-follower model. In the following, we set out the macroscopic control problem (without pointwise control constraints).

{Prob}

Let $\mathbf{g}$ be the state functions and let $u$ be the associated control function, we wish to solve the following optimisation problem:

[TABLE]

subject to

[TABLE]

with the initial conditions,

[TABLE]

and subject to no-flow boundary conditions.

A typical form for the cost functional $J$ is a quadratic cost,

[TABLE]

where $U$ is the set of admissible controls and $\beta>0$. Minimisation of this cost functional aims to concentrate opinions around a given desired state $w_{dF}$ while also minimising the energy used by the control variable to influence the state variable. Other choices of the cost functional are investigated in Section 6.1. In the following we consider the set of admissible controls

[TABLE]

for some $M>0$.

We can adapt results from [3] to establish the existence and uniqueness of solutions to (8) as well as the existence of optimal controls for Problem 4.1. We denote the dual space of $H^{1}(\mathcal{I})$ by $H^{-1}(\mathcal{I})$ and the duality pairing by $\langle\cdot,\cdot\rangle_{{H^{-1},H^{1}}}:H^{-1}(\mathcal{I})\times H^{1}(\mathcal{I})\to\mathbb{R},$ $\langle g,\varphi\rangle_{{H^{-1},H^{1}}}:=g(\varphi)$.

Definition 1**.**

For given $T>0$, functions $g_{L},g_{F}:\mathcal{I}\times(0,T)\to[0,\infty)$ are weak solutions of (8) if and only if

$g_{i}\in L^{2}(0,T;H^{1}(\mathcal{I}))$, $\tfrac{\partial g_{i}}{\partial t}\in L^{2}(0,T;H^{-1}(\mathcal{I}))$,  $i=L,F,$ 2. 2.

For any $\varphi\in L^{2}(0,T;H^{1}(\mathcal{I}))$,

[TABLE]

Theorem 1**.**

For given $T,M>0$, let $g_{L}^{0},g_{F}^{0}\in L^{2}(\mathcal{I})$ and $u\in U$. Then there exist unique weak solutions $g_{L},g_{F}$ to (8) in the sense of Definition 1.

Proof.

The existence and uniqueness for the leader equation (6) follows from Theorem 2.4 in [3]. With $g_{L}$ given, we can follow along the lines of the same proof for the follower equation (7), the additional term in (7) arising from the leader-follower interaction can be estimated similarly as the other interaction term, noting $g_{L}\in L^{2}(0,T;H^{1}(\mathcal{I})).$ ∎

Theorem 2**.**

For given $T,M>0$, let $\mathbf{g}^{0}=(g_{L}^{0},g_{F}^{0})\in L^{2}(\mathcal{I})$. Then there exist an optimal control $\bar{u}\in U$ with associated optimal state $\mathbf{\bar{g}}=(\bar{g}_{L},\bar{g}_{F})$.

Proof.

The proof uses standard arguments, following similar steps as in Theorem 2.5 in [3], presented here in concise form. Thanks to Theorem 1, for each $u\in U$ there exist unique solutions $g_{L},g_{F}$ to (8). Noting that the cost functional is bounded from below by zero, we can consider minimising sequences $\mathbf{g}^{j}=(g_{L}^{j},g_{F}^{j})$ and $u^{j}$ with

[TABLE]

We can employ Banach-Alaoglu’s theorem to extract converging subsequences,

[TABLE]

Noting that $U$ is weak-$\ast$ closed, we have $\bar{u}\in U.$ The above regularity is sufficient to pass to the limit in the weak formulation of (8) and to conclude that $\mathbf{\bar{g}}=(\bar{g}_{L},\bar{g}_{F})$ solves (8). Semicontinuity of the cost functional allows to conclude

[TABLE]

therefore $\bar{u}\in U$ is the optimal control with associated optimal state $\mathbf{\bar{g}}$. ∎

4.2 First-order optimality conditions

Using the formal Lagrangian method for Problem 4.1, by introducing the adjoint variable $\mathbf{p}=(p_{L},p_{F})^{\top}$ as our Lagrange multipliers we construct the formal Lagrangian,

[TABLE]

for which we must find the stationary point, $(\bar{\mathbf{g}},\bar{u},\bar{\mathbf{p}})$. To this end we compute the Fréchet derivatives of $\mathcal{L}$ with respect to $\mathbf{p}$, $\mathbf{g}$ and $u$, with the directions $\phi$, $\psi$ and $v$ respectively. That is we look for those functions $(\bar{\mathbf{g}},\bar{u},\bar{\mathbf{p}})$ such that,

[TABLE]

This leads to the weak form of the state system, $\nabla_{\mathbf{p}}\mathcal{L}\phi=0$,

[TABLE]

the weak form of the adjoint system, $\nabla_{\mathbf{g}}\mathcal{L}\psi=0$,

[TABLE]

and finally the optimality condition, $\nabla_{u}\mathcal{L}v=0$,

[TABLE]

where we omit the bars for clarity. We can now write (12)–(14) in their corresponding strong forms. These strong forms are somewhat unwieldy and therefore, we make a small assumption on the cost-functional $J$, to present the above system in more readable way – that is we assume that the cost functional will have the general form $J(\mathbf{g},\mathbf{u})=J_{T}[\mathbf{g}]+J_{s}[\mathbf{g},\mathbf{u}]$, that is, $J$ is separated into a time dependent part $J_{s}$ and a final time part $J_{T}$. For example the typical choice for $J$ stated in Section 3.1, the decomposition is $J_{T}[\mathbf{g}]=1/2\left\lVert\mathbf{g}(T)-\mathbf{g}_{\mathcal{I}}\right\rVert_{L^{2}(\mathcal{I};\real^{2})}^{2}$, $J_{s}[\mathbf{g},\mathbf{u}]=1/2\left\lVert\mathbf{u}\right\rVert_{U}^{2}$. We then write these equations in the strong form beginning with the state system,

[TABLE]

with initial condition $\mathbf{g}(0,w)=\mathbf{h}_{0}(w)$. We then write the strong form of the adjoint system,

[TABLE]

with the final time condition $\mathbf{p}(T,w)=\nabla_{\mathbf{g}}J_{T}[\mathbf{g}]$ and finally the strong form of the optimality condition,

[TABLE]

with the associated boundary conditions. In these equations, we required writing some extra terms, specifically,

[TABLE]

These terms were computed using a similar method to that in [3].

5 Numerical methods

In order to solve the optimality system from Section 4 numerically, we use a gradient-type sweeping algorithm [22]. This involves making an initial guess for the control function, $u_{-1}$, and then creating a sequence of functions $(\mathbf{g}_{i},u_{i},\mathbf{p}_{i})$ using the following algorithm: {Alg} With the initial guess $u_{-1}$ for the control function, set $n=0$,

1. [(1)]

Solve the state system $e(\mathbf{g_{n}},u_{n-1})=0$ forwards in time. 3. 2.

With the known $\mathbf{g}_{n}$, solve the adjoint system $e^{*}(\mathbf{g}_{n},u_{n-1},\mathbf{p}_{n})$ backwards in time. 4. 3.

Update the control function using $u_{n}=u_{n-1}-\nu\varepsilon(\mathbf{g}_{n},u_{n-1},\mathbf{p}_{n})$, for some chosen step size $\nu>0$. 5. 4.

Set $n=n+1$, return to step 1. until convergence of the cost functional $J$.

When solving the state system numerically some additional challenges present themselves. Since the solution of the state equation is a probability density function we require that $\mathbf{g}_{n}$ is unconditionally positive and the mass is conserved – $\int_{\mathcal{I}}g_{L,n}(t,w)\ \text{d}w=\int_{\mathcal{I}}g_{F,n}(t,w)\ \text{d}w=1$ – for all $s$. Additionally since the compromise term $A[\mathbf{g}_{n}]\mathbf{g}_{n}$ is non-linear and non-local special care must be taken when choosing a numerical method.

To combat these challenges, we use a Chang-Cooper method to discretise in space – which is second-order consistent in space and preserves quasi-steady states – and the modified Patankar-Runge-Kutta (MPRK) scheme to discretise in time – which is unconditionally positive [27]. One main advantage of this method over an explicit Euler method is the lack of constraints on the time step for convergence.

We partition our space axis with $L$ equidistant nodal points over $\mathcal{I}$, that is $w_{i}=-1+i\Delta w$ with $\Delta w=2/L$, and we partition our time axis with an equidistant mesh over $[0,T)$ with some $\Delta s$ chosen and $Q=\lfloor T/\Delta s\rfloor$ nodal points. We write a lower index to denote the spatial node value and an upper index to denote the temporal node value, that is $\mathbf{g}_{n}(w_{i},s^{t})=\mathbf{g}_{n,i}^{t}$. At each iteration in Algorithm 1 the values of $n$ are fixed, we omit these for ease of reading.

We now discuss the specifics of the discretisation in more detail. Since we are using a Chang-Cooper scheme for discretisation in space, we consider the semi-discretisation for (15),

[TABLE]

for the numerical fluxes,

[TABLE]

with

[TABLE]

With the semi-discretisation in space, we can write this equation in a production-destruction formulation, that is,

[TABLE]

with the corresponding initial condition and $\mathbf{g}$ representing the vector of time dependent functions and $p_{i,j}$ and $d_{i,j}$ representing the contribution to nodal point $i$ from nodal point $j$ – the production – and the loss from nodal point $j$ to nodal point $i$ – the destruction. Note that $p_{i,j}=d_{j,i}$. Since the solution of our state system are probability density functions we require that our scheme be conservative which allows us to narrow down how we can choose the functional $p$ and $d$, namely we wish for $p_{i,i}=d_{i,i}=0$ since this leads to

[TABLE]

which is equivalent to the conservation of mass for the semi-discretised functions $\mathbf{g}_{i}$.

We use the proposed values for $p$ and $d$ presented in [27] – where a more in depth description of this method can be found – namely,

[TABLE]

Finally, we can use the Patankar trick along with a second-order Runge-Kutta scheme for the time integration of (18),

[TABLE]

For the adjoint problem (16) the terms which require additional care are the non-local terms, $A[\mathbf{g}]\text{d}\mathbf{p}/\text{d}w$ and $A^{*}[\mathbf{g},\text{d}\mathbf{p}/\text{d}w]$. We use an upwind scheme to discretise the adjoint problem in space and discretise in time using a modified forward Euler scheme. The slight modification to the forward Euler scheme comes from treating the $A^{*}$ integral term in (16) explicitly with respect to $\mathbf{p}$. That is,

[TABLE]

which, after isolating the $\mathbf{p}^{t}$ terms on the left hand side and the other terms on the right hand side, yields,

[TABLE]

with

[TABLE]

Finally, the update of the control function is as described in Algorithm 1, evaluated at all space and time grid points.

6 Numerical experiments

In this section, we show the results of some numerical experiments, exploring the behaviour of these models using different compromise functions $P$ and cost functionals $J$. We take the initial guess of the control function in all experiments as $u\equiv 0$ for all $(w,t)\in\mathcal{I}\times[0,T]$, and define the state at time $t=0$ as localisation functions described in terms of its leader and follower components as,

[TABLE]

where $R_{L},R_{F}$ are the radii of the graph and $c_{L},c_{F}$ are the centres. In this equation, $\Lambda_{L}$ and $\Lambda_{F}$ are the constants which ensure,

[TABLE]

and $k$ is the sharpness factor of the hyperbolic tangent function. In particular we have chosen $R_{L}=R_{F}=0.85$, $c_{L}=c_{F}=0$, $k=10$.

For all the following experiments, we partition our opinion axis with $L=80$ equidistant nodal points in $\mathcal{I}$, that is $\Delta w=2/80=0.025$. Due to the absence of any time-stepping restriction in the MPRK method, we are free to choose a relatively large $\Delta s$ and we make the choice of $\Delta s=5\Delta w$ to reduce the time to run the simulations.

In Table 1, we report the parameters used in the numerical experiments for the Fokker-Planck modelling.

6.1 Investigating the choice of cost functional $J$

In this numerical experiment, we explore the behaviour of our model using some example cost functionals $J$ and compare these to the uncontrolled dynamics – a top view of the time evolution of which are presented in Figure 2. For this we use a bounded confidence model for our compromise functions, $P_{L}=P_{F}=\tilde{P}=P$ as follows,

[TABLE]

for $r=0.5$.

The first example cost functional is similar to that from [3, 8], which mirrors a common type of cost functional for the binary controlled dynamics (1)–(3),

[TABLE]

where $w_{dL},w_{dF}$ are some desired opinions that we wish the densities $g_{L}$ and $g_{F}$ to be centred around respectively.

We present the resulting states for the cost functional (22) in Figure 3 using $w_{dL}=w_{dF}=w_{d}=-0.5$. As we can see from these figures, the controlled states, mainly cluster around the desired position $w_{d}$ in this case, whereas the uncontrolled states form two distinct and equal clusters. In the follower dynamics we see a peak forming, centred at $w=0.8$, which – when looking at the time evolutions of the state – appears to be a concentration of followers that are outside the confidence radius $r=0.5$ of the leader species in the early time steps. This might be due to the fast convergence of the leaders to the desired state $w_{d}$, thus “leaving” some followers behind.

To explore this slightly more, we examine a variation of cost functional (22), where the cost functional depends only on the follower density rather than depending on both the leader and follower densities,

[TABLE]

From the results presented in Figure 4, we see the leader species follows a very similar strategy as in the previous experiment, however the follower species still forms a peak at around $w=0.8$, although this peak has less mass. A more illustrative example of this change in strategy for the leader species is where we modify the nature of the interactions between the leaders and followers. In this example we take

[TABLE]

known as a Sznajd-type interaction [24] for some constant $b$ and where $\tilde{P}$ is the compromise function that appears in the follower side of the leader-follower interaction (2).

Since this compromise function only depends on the interacting individuals’ opinions, this type of interaction can be seen as the propensity for individuals to change their opinion, where individuals with extreme opinions – close to $-1$ and $1$ – are influenced little by interactions and individuals with central opinions – close to [math] – are influenced greatly by the interactions. The constant $b$ is used to describe the weight of the interaction and in general we allow $b\in[-1,1]$ to ensure that $|\tilde{P}|<1$, with $b>0$ representing a concentration of the follower density around the leader density and $b<0$ representing the separation of the followers from the leaders.

In this numerical experiment we use $b=-1$, to ensure that either the leader species or the follower species can converge on the desired opinion but not both. Recall that the interactions within the leader and follower species are still using the bounded confidence-type interactions – interactions only occur when the opinions of individuals are close together. This leads to radically different strategies for the leader species depending on the choice of cost functional.

The results presented in Figure 5 use the cost functional (22) and we can see that although the leader species concentrates around $w_{d}=-0.5$ the follower species does not vary from the uncontrolled dynamics by much. In these dynamics, we see an increase in the early values for $J$ – which is due to the uncontrolled nature of the follower species. The density of the leader species starts to approach $w_{d}$ at around sweeping iteration 6, however the follower density is still at both the extremes $+1$ and $-1$. As the leaders move from the positive half of the domain to the negative half, the followers that start with an opinion around $w=0$ move to the boundary at $-1$ rather than that at $+1$. This then leads to a new leader strategy to minimise in this new regime. This effect can be seen in the presented time evolution for the leader density, Figure 5(c), where in earlier time steps, a concentration of the leaders close to $+1$ causes the follower density to form a larger peak at the $-1$ boundary.

When using the cost functional (23) – presented in Figure 6 – the leader strategy changes dramatically to cause the concentration of the followers about the desired opinion.

Finally, we look at the cost functional

[TABLE]

where we return to using the bounded confidence model, that is $P_{L}=P_{F}=\tilde{P}=P$ as in the first experiment. In the experiment presented here, the target function is defined with \mathbf{g}_{\mathcal{I}}(w)=\big{(}g_{\mathcal{I}}(w),g_{\mathcal{I}}(w)\big{)} where $g_{\mathcal{I}}$ is constructed by setting specific values for the $w\in\{-0.5,0,0.5\}$ and then using a Lagrangian polynomial interpolation of these points with a repeated root at $w=-1$ and $w=1$, i.e. $g_{\mathcal{I}}(-1)=g_{\mathcal{I}}(1)=0$. The specific values for this experiment were chosen as $g_{\mathcal{I}}(-0.5)=0.5$, $g_{\mathcal{I}}(0)=0.25$ and $g_{\mathcal{I}}(0.5)=1$.

As can be seen from Figure 7, the leader species changes dramatically from the uncontrolled dynamics to closer match the target density $g_{\mathcal{I}}$. Although a less drastic change has happened to the follower density, a clear asymmetry has appeared. Since the follower dynamics only see the control through their interactions with the leader species, we expect that the leader species would minimise the controlled dynamic, then the follower dynamic follows which reduces the total cost. In order to see a closer fit of the follower dynamics to the target, we again modify (24) to use only the follower density,

[TABLE]

Once again – as presented in Figure 8 – we see a large change in the leader density but a closer match of the follower dynamics to the desired final time density.

Declarations

• •

Funding: This work was partially supported by a Royal Society International Exchanges grant (Ref. IES/R3/213113). O.W. is supported by the Warwick Mathematics Institute Centre for Doctoral Training in Mathematics, and gratefully acknowledges funding from the University of Warwick.

• •

Conflict of interest: On behalf of all authors, the corresponding author states that there is no conflict of interest.

Appendix A Derivation of controlled Boltzmann system

We consider the evolution of the random variables that describe the opinion of an individual in the leader species $X_{L}$, in a the small time interval $[t,t+\Delta t]$ for some small $\Delta t$. We assume this time step is chosen such that $X_{L}$ interacts at most once, and the interactions that can occur will be of type (1), (2) or $X_{L}$ does not interact in this time interval. We can now introduce the random variable, $T$, with probability mass function,

[TABLE]

with $\alpha,\beta>0$ weighting the probability of interaction with the leader and follower respectively. Here $T=0$ represents the interaction (1), $T=1$ represents (2) and $T=2$ represents no interaction. Labelling the other potential interacting individuals as $Y_{L}$ and $Y_{F}$ as the leader and follower respectively, we can write the opinion $X_{L}$ at time $t+\Delta t$ as,

[TABLE]

with $X^{*}_{L}[X_{L},Y_{L}]$ representing the position of $X_{L}$ under the interaction (1), and $X_{L}^{*}[X_{L},Y_{F}]$ the position of $X_{L}$ under the interaction (2). Note that in our case $X_{L}^{*}[X_{L},Y_{F}]=X_{L}$.

Now considering $\left\langle\phi(X_{L}(t+\Delta t))\right\rangle$ – where $\left\langle\cdot\right\rangle$ is the mean with respect to $X_{L}$, $Y_{L}$, $Y_{F}$ and $\eta_{L}$ and $\phi\in C_{c}^{2,\delta}(\mathcal{I};\real)$ – and using (26), we obtain,

[TABLE]

through straightforward calculation, we get,

[TABLE]

we now apply the limit as $\Delta t\rightarrow 0$,

[TABLE]

Note here that the individual $Y_{L}$ is also a member of the leader species and thus the density which describes the distribution of $X_{L}$ also describes that of $Y_{L}$. Writing the expectations of $X_{L}$, $Y_{L}$ and $Y_{F}$ explicitly yields the weak form of the Boltzmann-type equation for the leader species,

[TABLE]

Notice that the second term on the right hand side of (29) will cancel out, since the post-interaction opinion of the leader with a follower results in no change to the leader opinion, that is $\phi(w^{*})-\phi(w)=0$. We can then write the strong equation,

[TABLE]

defining the $Q_{L}$ operator as the strong form of the remaining term on right hand side of (29).

Using a slight modification of the above method for the follower species, we derive the remaining equations of the system,

[TABLE]

Once again, since $X_{F}$ and $Y_{F}$ are both members of the follower species, the densities which describe their distributions are the same, writing the expectation in (31) explicitly, we obtain the weak form of the Boltzmann-type equation for the follower species,

[TABLE]

with the strong form of this equation,

[TABLE]

The interaction operators $Q_{L}$, $Q_{FL}$, and $Q_{F}$ in (30) and (33), are expressed in weak form as,

[TABLE]

where $\phi$ are test functions from $C_{c}^{2,\delta}(\mathcal{I};\real)$, that are continuous and compactly supported functions on $\mathcal{I}$. In (30) and (33), the values of $\alpha$, $a$, $b$ denote the interaction frequencies, which we can write as relaxation times of the Boltzmann-type equations using $\alpha=1/\tau_{LL}$, $a=1/\tau_{FF}$ and $b=1/\tau_{FL}$.

Appendix B Derivation of controlled Fokker-Planck system

We derive the coupled system of non-linear, non-local Fokker-Planck equations for this model, starting with the leader equation by considering the first order Taylor’s expansion of $\phi(w^{*})$ about $w$,

[TABLE]

with $\zeta$ and $R$ defined as,

[TABLE]

where $\tilde{w}=\kappa w+(1-\kappa)w^{*}$ for $\kappa\in[0,1]$. Applying the same Taylor’s expansion for $\phi(v^{*})$ about $v$, substituting into $Q_{L}$ and then iterating the integral on the right hand side yields,

[TABLE]

where,

[TABLE]

Now considering the time-scaling $s=\gamma_{L}t$ and the transformed density functions $g_{L}(w,s)=f_{L}(w,t)$ and taking the quasi-invariant limit of (30) – that is $\sigma,\gamma_{L}\rightarrow 0$ while keeping the value $\sigma^{2}/\gamma_{L}$ constant throughout. Notice that under the time scaling and quasi-invariant limit, $R$ goes to zero (see [1]) and in the $\zeta$ term, we still have a factor of $\gamma_{L}$ even after dividing through by $\gamma_{L}$ thus this term also goes to zero. We then write the strong form of the Fokker-Planck-type equation for the leader species,

[TABLE]

with $\lambda_{L}=M_{L}\sigma^{2}/\gamma_{L}$ and $M_{L}=\int_{\mathcal{I}}f_{L}(w,t)\,dw$, subject to no-flow boundary conditions.

Applying the same process for (33), with the scaling $\gamma_{L}=\alpha_{LF}\gamma_{F}$ and letting $\lambda_{F}=M_{F}\sigma^{2}/\gamma_{F}$ with $M_{F}=\int_{\mathcal{I}}f_{F}(w,t)\,dw$, yields the follower Fokker-Planck-type equation,

[TABLE]

subject to no-flow boundary conditions, where:

[TABLE]