# Contrastive Brain Network Learning via Hierarchical Signed Graph Pooling
Model

Haoteng Tang, Guixiang Ma, Lei Guo,
Xiyao Fu, Heng Huang, Liang Zhan† H. Tang, L. Guo, X. Fu, H. Huang and L. Zhan
are with the Department of Electrical and Computer Engineering, University of
Pittsburgh, Pittsburgh, PA, 15260, USA (e-mail: {haoteng.tang, Lei.guo, xy_fu,
heng.huang, liang.zhan}@pitt.edu)G. Ma is with the Department of Computer
Science, University of Illinois at Chicago, Chicago, IL 60607, USA
(e-mail:guixiang.ma@intel.com). This work was done before Dr. Ma joined
Intel.Liang Zhan is the corresponding author (denoted by ${\dagger}$)

###### Abstract

Recently brain networks have been widely adopted to study brain dynamics,
brain development and brain diseases. Graph representation learning techniques
on brain functional networks can facilitate the discovery of novel biomarkers
for clinical phenotypes and neurodegenerative diseases. However, current graph
learning techniques have several issues on brain network mining. Firstly, most
current graph learning models are designed for unsigned graph, which hinders
the analysis of many signed network data (e.g., brain functional networks).
Meanwhile, the insufficiency of brain network data limits the model
performance on clinical phenotypes predictions. Moreover, few of current graph
learning model is interpretable, which may not be capable to provide
biological insights for model outcomes. Here, we propose an interpretable
hierarchical signed graph representation learning model to extract graph-level
representations from brain functional networks, which can be used for
different prediction tasks. In order to further improve the model performance,
we also propose a new strategy to augment functional brain network data for
contrastive learning. We evaluate this framework on different classification
and regression tasks using the data from HCP and OASIS. Our results from
extensive experiments demonstrate the superiority of the proposed model
compared to several state-of-the-art techniques. Additionally, we use graph
saliency maps, derived from these prediction tasks, to demonstrate detection
and interpretation of phenotypic biomarkers.

###### Index Terms:

Signed Graph Learning, Hierarchical Graph Pooling, Contrastive Learning, Brain
Functional Networks, Data Augmentation, Interpretability.

##  I Introduction

Understanding brain organizations and their relationship to phenotypes (e.g.,
clinical outcomes, behavior or demographical variables, etc.) are of prime
importance in the modern neuroscience field. One of important research
directions is to use non-invasive neuroimaging data (e.g., functional magnetic
resonance imaging or fMRI) to identify potential imaging biomarkers for
clinical purposes. Most previous research focuses on voxel-wise and region-of-
interests (ROIs) imaging features [1, 2, 3]. However, evidences show that most
of these clinical or behavior phenotypes are the outcomes of interactions
among different brain regions. Therefore, brain networks attract more and more
attention for the purpose of phenotype predictions [4, 5, 6]. Additionally,
compared to traditional neuroimaging features, brain network has more
potential to gain interpretable and system-level insights into phenotype-
induced brain dynamics [7]. A brain network is a 3D brain graph model, where
graph nodes represent the attributes of brain regions and graph edges
represent the connections (or interactions) among these regions.

Many studies have been conducted to analyze brain networks based on the graph
theory, however, most of these studies focus on pre-defined network features,
such as clustering coefficient, small-worldness [8, 9, 10, 11, 12]. This may
be sub-optimal since these pre-defined network features may not be able to
capture the characteristics of the whole brain network. However, the whole
brain network is difficult to be analyzed due to the high dimensionality. To
tackle this issue, Graph Neural Network (GNN), as one of embedding techniques,
has gained increasing attentions to explore biological characteristics of
brain network-phenotype associations in recent years [13, 14, 15]. GNN is a
class of deep neural networks that can embed the high-dimensional graph
topological structures with graph node features into low dimensional latent
space based on the information passing mechanism [16, 17, 18]. A few studies
proposed different GNNs to embed the nodes in brain networks and applied a
global readout operation (e.g., global mean or sum) to summarize all latent
node features as the whole brain network representation for downstream tasks
(e.g., behavior score regression, clinical disease classification) [14, 15,
19]. However, the message passing of GNNs is inherently ‘flat’ which only
propagates information across graph edges and is unable to capture
hierarchical structures rooted in graphs which are crucial in brain functional
organizations [20, 21, 22, 23]. To address this issue, many recent studies
introduce hierarchical GNNs, including node embedding and hierarchical graph
pooling strategies, to embed the whole brain network in a hierarchical manner
[20, 24, 25, 26, 27].

Although GNNs have achieved great progresses on brain network mining, three
issues should be addressed:

  * •

Most current GNNs are designed for unsigned graphs in which all graph nodes
are connected via non-negative edges (i.e., edge weights are in the range of
$[0,\infty)$). However, signed graphs are very common in brain research (e.g.,
functional MRI-derived brain networks or brain functional networks).
Therefore, signed graph embedding models are valuable.

  * •

Brain network data, compared with other types of network data, is insufficient
since the data collection is very expensive and time consuming. This may limit
the model performance on prediction tasks in a way.

  * •

Most current GNNs on brain network studies are not interpretable, and thus are
incapable to provide biological explanations or heuristic insights for model
outcomes. This is mainly due to the black-box nature of the neural networks.

To tackle the first issue, a few recent studies proposed signed graph
embedding models based on the balance-theory [28, 29, 30, 31]. The balance-
theory, motivated by human attitudes in social networks, is used to describe
the node relationship in signed graphs, where nodes connected by positive
edges are considered as ‘friends’ otherwise are considered as ‘opponents’.
Meanwhile, the balance-theory also defines $4$ higher-order relationships
among graph nodes: (1) the ‘friend’ of ‘friend’ is ‘friend’, (2) the
‘opponent’ of ‘friend’ is ‘opponent’, (3) the ‘friend’ of ‘opponent’ is
‘opponent’ and (4) the ‘opponent’ of ‘opponent’ is ‘friend’. These definitions
are accorded with nodal relationships in the functional brain network, which
indicates that the balance theory might be applicable in functional brain
network embedding. However, existing signed graph embedding models focus on
embedding graph nodes with signed edges into latent features without
considering the hierarchical structures in graphs, which may not facilitate
the whole graph representation learning and the graph-level tasks (i.e.,
clinical disease classification based on whole brain networks). To address
this issue, we propose a hierarchical graph pooling module on signed graphs
based on the information theory and extend the current methods to a
hierarchical signed graph embedding model.

To address the second issue, we propose a data augmentation strategy to
augment functional brain networks. Meanwhile, we introduce the graph
contrastive learning architecture, where contrastive graph samples are
generated by the proposed augmentation strategy, to boost the model
performance on prediction tasks. The data augmentation aims at creating
reasonable data samples, by applying certain transformations, which are
similar to the original ones. For example, image rotation and cropping are
common transformations to generate new samples in image classification tasks
[32, 33, 34]. In graph structural data, a few studies proposed to utilize
graph perturbations (i.e., add/drop graph nodes, manipulate graph edges) and
graph view augmentation (e.g., graph diffusion) to generate contrastive graph
samples from different views [35, 36, 37, 38]. These strategies, although
boosting the model performance on large-scale benchmark datasets (e.g., CORA,
CITESEER, etc.), may not be suitable to generate contrastive brain network
samples. On the one hand, each node in brain networks represents a defined
brain region with specific brain activity information so that the brain node
can not be arbitrarily removed or added. On the other hand, add/drop
operations on brain network may lead to unexpected model outcomes which are
difficult to explain and understand from biological views. Therefore, we
generate the augmented brain functional networks directly from fMRI BOLD
signals, where the generated samples are similar and the biological structure
is maintained.

As for the last issue, our proposed graph pooling module is interpretable by
nature. Previous studies indicated that brain networks are hierarchically
organized by some regions as neuro-information hubs and peripheral regions,
respectively [39, 40, 41]. Within our graph pooling module, an information
score is designed to measure the information gain for each brain node and only
top-$K$ nodes with high information gains will be preserved as brain
information hubs while the information of other peripheral brain nodes will be
aggregated onto these hubs. Hence, the proposed pooling module can be
interpreted as a brain information hubs generator. Apparently, the outcome of
this pooling module is a subgraph of the original brain network without any
new nodes. Therefore, yielded subgraph nodes can be regarded as potential
biomarkers to provide heuristic biological explanations for tasks.

Our main contributions are summarized as follow:

  * •

We propose a hierarchical signed graph representation learning (HSGRL) model
to embed brain functional networks and we apply the proposed model on multiple
phenotype prediction tasks.

  * •

We propose an augmentation strategy for fMRI-derived brain network data. To
further boost the model performance, we build up a contrastive learning
framework with the proposed HSGPL model, where the contrastive samples are
generated by the designed augmentation strategy.

  * •

The proposed HSGPL model is interpretable which yields heuristic biological
explanations.

  * •

Extensive experiments are conducted to demonstrate the superiority of our
method. Moreover, we draw graph saliency maps for clinical tasks, to enable
interpretable detection of phenotype biomarkers.

##  II Related Works

###  II-A Graph Neural Networks and Brain Network Embedding

GNNs are generalized deep learning architectures which are broadly utilized
for graph representation learning in many fields (e.g., social network mining
[42, 43], molecule studies [44, 45] and brain network analysis [46]). Most
existing GNN models (e.g., GCN [16], GAT [17], GraphSage [47]) focus on node-
level representation learning and only propagate information across edges of
the graph in a flat way. When deploying these models on graph-level tasks
(e.g., graph classification, graph similarity learning, [48, 49, 50, 51]), the
whole graph representations are obtained by a naive global readout operation
(e.g., sum or average all node feature vectors). However, this may lead to
poor performance and low efficiency in graph-level tasks since the
hierarchical structure, an important property that existed in graphs, is
ignored in these models. To explore and capture hierarchical structures in
graphs, a few hierarchical graph pooling strategies are proposed to learn
representations for the whole graph in a hierarchical manner [24, 20, 52, 53,
25]. Traditional methods to extract brain network patterns are based on graph
theory [8, 9, 10, 11, 12] or geometric network optimization [54, 55, 56, 57].
A few recent studies [14, 15, 58] introduce GNNs to discover brain patterns
for phenotypes predictions. However, hierarchical structures in brain networks
are not considered in these models, which limits the model performance in a
way. Recently, a few hierarchical brain network embedding models are
proposed[59, 26].

However, all the aforementioned GNNs are designed for unsigned graph
representation learning. A few recent studies are proposed to handle the
signed graphs, however, they only consider the node-level representation
learning [29, 31, 60, 61]. In this work, we design a signed graph hierarchical
pooling strategy to extract graph-level representations from brain functional
networks.

###  II-B Interpretable Graph Learning Model

Generally, the mechanism about how GNNs embed the graph nodes can be explained
as a message passing process, which includes message aggregations from
neighbor nodes and message (non-linear) transformations [18, 62, 26]. However,
most current hierarchical pooling strategies are not interpretable [20, 25,
24]. A few recent studies try to propose interpretable graph pooling
strategies to make the pooling module intelligible to the model users. Most of
these pooling strategies down-sample graphs relying on network communities
which are one of the important hierarchical structures that can be interpreted
[27, 26, 63]. For example, [26] proposed a hierarchical graph pooling neural
network relying on brain network community to yield interpretable biomarkers.
The hierarchical pooling strategy proposed in this work relies on the network
information hub which is another important hierarchical structure in brain
networks.

###  II-C Data Augmentation for Graph Contrastive Learning

Most current graph contrastive learning methods augment graph contrastive
samples by manipulating graph topological structures. For example, [36, 37]
generate the contrastive graph samples by dropping nodes and perturbing edges.
Other studies generate contrastive samples by changing the graph local
receptive field, which is named as the graph view augmentation [64, 35]. In
this work, we introduce the graph contrastive learning into brain functional
network analysis and generate contrastive samples from the fMRI BOLD signals.

##  III Preliminaries of Brain Functional Networks

We denote a brain functional network with $N$ nodes as $G=\\{V,E\\}=(A,H)$.
$V$ is the graph node set where each node (i.e., $v_{i},i=1,...,N$) represents
a brain region. $E$ is the graph edge set where each edge (i.e., $e_{i,j}$)
describes the connection between node $v_{i}$ and $v{j}$.
$A\in\mathbb{R}^{N\times N}$ is the graph adjacency matrix where each element,
$a_{i,j}\in A$, is the weight of edge $e_{i,j}$. $H\in\mathbb{R}^{N\times C}$
is the node feature matrix where $H_{i}\in H$ is the $i-th$ row of $H$
representing the feature vector of $v_{i}$. Let $B\in\mathbb{R}^{N\times D}$
be the fMRI BOLD signal matrix, where $D$ is the signal length. Generally, the
edge weight in the brain functional network can be computed from the fMRI BOLD
signal by $a_{i,j}=corr(b_{i},b_{j})$, where $b_{i}$ is the $i-th$ row of $B$
representing the BOLD signal of $v_{i}$ and $corr(\cdot)$ is the correlation
coefficient operator. Note that $a_{i,j}$ can be either positive or negative
value so that brain functional network is a signed graph. For each subject, we
use $\hat{}$ and $\check{}$ to denote a functional brain network contrastive
sample pair (i.e.,
$[\hat{G}=(\hat{A},\hat{H}),\check{G}=(\check{A},\check{H})]$).

##  IV Methodology

Figure 1: Diagram of the proposed contrastive graph learning framework (in the
bottom black box) with hierarchical signed graph representation learning model
(in the top black box) for functional brain network embedding and downstream
tasks (i.e., phenotype classification or regression).

In this section, we first propose a data augmentation strategy to generate
contrastive samples for brain functional networks. Secondly, we introduce our
proposed hierarchical signed graph representation learning (HSGRL) model with
node embedding and hierarchical graph pooling modules. Finally, we deploy the
contrastive learning framework on our proposed HSGRL model to yield the
representations for the whole graph, which can be applied to downstream
prediction tasks.

###  IV-A Contrastive Samples of Brain Functional Networks

The generation of contrastive samples aims at creating reasonable and similar
functional brain network pairs by applying certain transformations. Here we
propose a new strategy to generate the brain functional network contrastive
samples from fMRI BOLD signals. For each node $v_{i}$, we generate two sub-
BOLD-signals ($\hat{b}_{i}$ and $\check{b}_{i}$) by manipulating its original
bold signal $b_{i}$. Specifically, we use a window ($size=d$) to clamp the
$b_{i}$ from the signal head and tail, respectively:

$\displaystyle\hat{b}_{i}$ $\displaystyle=$ $\displaystyle b_{i}[d+1,\quad
d+2,\quad...,\quad D]$ $\displaystyle\check{b}_{i}$ $\displaystyle=$
$\displaystyle b_{i}[1,\quad 2,\quad...,\quad D-d]$ (1)

Obviously, $b_{i}\in\mathbb{R}^{1\times D}$, $\hat{b}_{i}$ and
$\check{b}_{i}\in\mathbb{R}^{1\times(D-d)}$. To keep the similarity between
$\hat{G}$ and $\check{G}$, we set the window size $d\ll D$. After we generate
a pair of sub-bold-signals, we can compute edge weights of the pairwise
contrastive brain functional network samples by:

$\displaystyle\hat{a}_{i,j}$ $\displaystyle=$ $\displaystyle
corr(\hat{b}_{i},\hat{b}_{j})$ $\displaystyle\check{a}_{i,j}$ $\displaystyle=$
$\displaystyle corr(\check{b}_{i},\check{b}_{j}),$ (2)

where $\hat{a}_{i,j}\in\hat{A}$ and $\check{a}_{i,j}\in\check{A}$ are the
weights of $e_{i,j}$ in two contrastive samples. We do not consider the
contrastive node features in this work, therefore $\hat{X}=\check{X}=X$. The
generated contrastive sample pairs are similar with same node features and
slightly different edge weights. We will show this similarity in section V-C.

###  IV-B Hierarchical Signed Graph Representation Learning Model

We present our Hierarchical Signed Graph Representation Learning (HSGRL) model
in Fig. 1. The HSGRL model includes Balanced and Unbalanced Embedding (BUE)
module and Hierarchical Graph Pooling (HGP) module.

####  IV-B1 BUE module

The balance theory is broadly used to analyze the node relationships in signed
graphs. The theory states that given a node $v_{i}$ in a signed graph, any
other node (i.e., $v_{j}$) can be assigned into either balanced node set or
unbalanced node set to $v_{i}$ regarding to a path between $v_{i}$ and
$v_{j}$. Specifically, if the number of negative edges are even in the path
between $v_{i}$ and $v_{j}$, then $v_{j}$ belongs to the balanced set of
$v_{i}$. Otherwise, $v_{j}$ belongs to the unbalanced set of $v_{i}$. The
balance theory indicates that:

  * •

Each graph node, $v_{j}$, can belong to either the balanced or unbalanced node
set of a given target node $v_{i}$.

  * •

The path between $v_{i}$ and $v_{j}$ determines the balance attribute of
$v_{j}$.

Motivated by this, we adopt the idea of signed graph attention networks from
[31] to embed brain functional network nodes to generate latent node features
with balanced and unbalanced components:

$\displaystyle X^{B},X^{U}=F_{sign}(A,H)$ (3)

where $F_{sign}(\cdot)$ is the signed graph attention encoder [31]. $X^{B}$
and $X^{U}$ are the node balanced and unbalanced components of node latent
features, respectively. We fuse the two feature components as the node latent
features by:

$\displaystyle X=[X^{B}\|X^{U}],$ (4)

where $[||]$ denotes concatenate operation.

####  IV-B2 Hierarchical Signed Graph Pooling

As shown in Fig 1, the proposed Hierarchical Graph Pooling (HGP) module
consists of $4$ steps including: (A) information scores computation, (B) Top-K
informative hubs selection, (C) features aggregation and (D) graph pooling.

Information Score Computation: The information score of each node is also
considered to contain balanced and unbalanced components to measure the
information quantity that each node gains from balanced node set and
unbalanced node set, respectively. We first split the signed graph (i.e., with
adjacency matrix as $A$) into positive sub-graph (with adjacency matrix as
$A_{+}$) and negative one (with adjacency matrix as $A_{-}$). Then we utilize
Laplace normalization to normalize these two adjacency matrices as:

$\displaystyle\bar{A}_{+}$ $\displaystyle=$ $\displaystyle
D^{-\frac{1}{2}}_{+}A_{+}D^{-\frac{1}{2}}_{+}$ $\displaystyle\bar{A}_{-}$
$\displaystyle=$ $\displaystyle D^{-\frac{1}{2}}_{-}\lvert A_{-}\rvert
D^{-\frac{1}{2}}_{-},$ (5)

where $\bar{A}$ is the normalized adjacency matrix. $D_{+}$ and $D_{-}$ are
degree matrices of $A_{+}$ and $|A_{-}|$, respectively. Note that the i-th
line in $\bar{A}$, denoted by $\bar{A}_{i}$, represents the connectivity
probability distribution between $v_{i}$ and any other nodes. For each node
(i.e., $v_{i}$), we respectively define the balanced and unbalanced components
of information score (IS) by:

$\displaystyle IS^{B}_{i}$ $\displaystyle=$
$\displaystyle\left\|\bar{A}_{+,i:}^{\top}\otimes
X^{B}\right\|_{\tilde{L}_{1}}+\left\|\bar{A}_{-,i:}^{\top}\otimes
X^{U}\right\|_{\tilde{L}_{1}}$ $\displaystyle IS^{U}_{i}$ $\displaystyle=$
$\displaystyle\left\|\bar{A}_{+,i:}^{\top}\otimes
X^{U}\right\|_{\tilde{L}_{1}}+\left\|\bar{A}_{-,i:}^{\top}\otimes
X^{B}\right\|_{\tilde{L}_{1}},$ (6)

where $\left\|\cdot\right\|_{\tilde{L}_{1}}$ is line-wise $L_{1}$ norm, and
$\otimes$ is the scalar-multiplication between each line of two matrices.
$\top$ represents transpose of vector. Then the IS of $v_{i}$ can be obtained
by:

$\displaystyle IS_{i}=IS^{B}_{i}+IS^{U}_{i}.$ (7)

Top-$K$ Node Selection and Feature Aggregation: After we obtain the
information score for each brain node, we rank the IS and select $K$ brain
nodes, with top-$K$ IS values, as informative network hubs. For the other
nodes, we aggregate their features on the selected $K$ network hubs based on
the feature attention. Particularly, the feature attention between $v_{i}$ and
$v_{j}$ is computed by: $x_{i}x_{j}^{\top}$. We weighted add (i.e., set
feature attentions as weights) the feature of each unselected node to one of
hub features, where the attention value between these two nodes is the
biggest.

Graph Pooling After the feature aggregation, we down-scale the graph node by
removing all unselected nodes. In another word, only the selected top-$K$
network hubs as well as the edges among them will be preserved after graph
pooling. Since the functional brain network is a fully connected graph so that
no isolated node is existed in the down-scaled graph.

###  IV-C Contrastive Learning Framework with BUE and HGP

The contrastive learning framework with HSGRL is presented in Fig. 1. Assume
that we forward a pair of contrastive graph samples into the proposed HSGRL
model, we will obtain two node latent features, $\hat{X}$ and $\check{X}$
after the last pooling module. We first generate the graph-level
representations of two functional brain networks based on the latent node
features by a readout operator:

$\displaystyle\hat{X}_{G}=\sum_{i=1}^{N^{\prime}}\hat{x}_{i},\quad\check{X}_{G}=\sum_{i=1}^{N^{\prime}}\check{x}_{i},$
(8)

where $\hat{x}_{i}$ and $\check{x}_{i}$ are $i-th$ row of $\hat{X}$ and
$\check{X}$. $N^{\prime}(<N)$ is the number of nodes in the down-scaled graph
generated by the last pooling module.

####  IV-C1 Contrastive Loss

The normalized temperature-scaled cross entropy loss [65, 66, 67] is utilized
to construct the contrastive loss. In the framework training stage, we
randomly sample $M$ pairs from the generated contrastive graph samples as a
mini-batch and forward them to the proposed HSGRL model to generate
contrastive graph representation pairs (i.e., $\hat{X}_{G}$ and
$\check{X}_{G}$). We use $m\in\\{1,...,M\\}$ to denote the ID of the sample
pair. The contrastive loss of the $m-th$ sample pair is fomulated as:

$\displaystyle\ell_{m}=-log\frac{exp(\Phi(\hat{X}_{G}^{m},\check{X}_{G}^{m})/\alpha)}{\sum_{t=1,t\neq
m}^{M}exp(\Phi(\hat{X}_{G}^{m},\check{X}_{G}^{t})/\alpha)},$ (9)

where $\alpha$ is the temperature parameter. $\Phi(\cdot)$ denotes a
similarity function that:

$\Phi(\hat{X}_{G}^{m},\check{X}_{G}^{m})=\hat{X}_{G}^{m\top}\check{X}_{G}^{m}/\|\hat{X}_{G}^{m}\|\|\check{X}_{G}^{m}\|.$
(10)

The batch contrastive loss can be computed by:

$\mathcal{L}_{contrastive}=\frac{1}{M}\sum_{m=1}^{M}\ell_{m}$ (11)

####  IV-C2 Downstream Task and Loss Functions

We use an MLP to generate the framework prediction for both classification and
regression tasks. Specifically, the prediction can be generate by
$Y_{pred}=MLP([\hat{X}_{G}\|\check{X}_{G}])$. We use $NLLLoss$ and $L_{1}Loss$
as supervised loss functions ($\mathcal{L}_{supervised}$) of classification
and regression tasks, respectively. The whole framework can be trained in an
end-to-end manner by optimizing:

$\mathcal{L}=\mu_{1}\mathcal{L}_{supervised}+\mu_{2}\mathcal{L}_{contrastive},$
(12)

where $\mu_{1}$ and $\mu_{2}$ are the loss weights.

Figure 2: Visualization of the averaged adjacency matrices for original and
contrastive samples. The averaged contrastive sample pair is generated by
using a window size $d=10$.

##  V Experiments

###  V-A Datasets and Data Preprocessing

Two publicly available datasets were used to evaluate our framework. The first
includes $1206$ young healthy subjects (mean age $28.19\pm 7.15$, $657$ women)
from the Human Connectome Project (HCP) [68]. The second includes 1326
subjects (mean age $=70.42\pm 8.95$, $738$ women) from the Open Access Series
of Imaging Studies (OASIS) dataset [69]. Details of each dataset can be found
on their official websites 111https://www.oasis-brains.org
222https://wiki.humanconnectome.org. CONN [70] were used to preprocess fMRI
data and the preprocessing pipeline follows our previous publications [71,
72]. For HCP data, each subject’s network has a dimension of $82\times 82$
based on 82 ROIs defined using FreeSurfer (V6.0) [73]. For OASIS data, each
subject’s network has a dimension of $132\times 132$ based on the Harvard-
Oxford Atlas and AAL Atlas. We deliberately chose different network
resolutions for HCP and OASIS to evaluate whether the performance of our new
framework is affected by the network dimension or atlas.

###  V-B Implementation Details

We randomly split the entire functional brain network dataset into $5$
disjoint subsets for $5$-fold cross-validations in our experiments. The values
in the adjacency matrices ($\hat{A}$ and $\check{A}$) of brain functional
networks are within range of $[-1,1]$. We compute the kurtosis and skewness
values of the fMRI BOLD signals as the node feature matrices ($H$). We use the
Adam optimizer [74] to optimize the loss functions in our model with a batch
size of $128$. The initial learning rate is $1e^{-4}$ and decayed by
$(1-\frac{current\\_epoch}{max\\_epoch})^{0.9}$. We also regularized the
training with an $L_{2}$ weight decay of $1e^{-5}$. We set the maximum number
of training epochs as $1000$ and, following the strategy in [24],[75], stop
training if the validation loss does not decrease for $50$ epochs. The
experiments were deployed on one NVIDIA RTX A6000 GPU.

###  V-C Similarities of Contrastive Samples

We utilize the $L_{2}$ distance and Cosine Similarity to measure the
similarities of the adjacency matrices of contrastive brain networks. Here, we
set the window size $d=10$ to generate the contrastive adjacency matrices. The
inner-pair similarity is computed by
$\frac{1}{M}\sum_{m=1}^{M}\Psi(\hat{A}^{m},\check{A}^{m})$, and the inter-pair
similarity is computed by
$\frac{1}{M^{2}}\sum_{m=1}^{M}\sum_{t=1}^{M}\Psi(\hat{A}^{m},\check{A}^{t})$,
where $\Psi(\cdot)$ is the similarity function (i.e., $L_{2}$ distance or
Cosine Similarity). The inner-pair $L_{1}$ distances on HCP and OASIS data are
$0.1301$ and $0.0915$, respectively. The inner-pair Cosine Similarities on HCP
and OASIS data are $0.9283$ and $0.9466$, respectively. The inter-pair $L_{1}$
distances on HCP and OASIS data are $0.2925$ and $0.3137$, respectively. The
inter-pair Cosine Similarities on HCP and OASIS data are $0.7311$ and
$0.7014$, respectively. We visualize the averaged adjacency matrics on HCP
data in Fig. 2 to show their similarities. The original sample is generated by
using the whole fMRI BOLD signal (i.e., $d=0$).

###  V-D Classification Tasks

Table I: Classification accuracy with s.t.d values under 5-fold cross-validation on gender classification, zygosity classification and AD classification tasks. The values in bold show the best results. Method | HCP | OASIS
---|---|---
Gender | Zygosity | AD
Acc. | Pre. | F1. | Acc. | Macro-F1. | Acc. | Pre. | F1.
t-BNE | 63.84(2.09) | 64.17(1.90) | 63.264(2.12) | 37.19(2.65) | 39.67(3.04) | 61.26(2.31) | 63.58(2.06) | 62.05(1.97)
mCCA-ICA | 61.21(4.03) | 63.11(3.75) | 62.20(3.59) | 35.51(4.64) | 38.71(3.34) | 63.37(1.98) | 62.06(2.12) | 64.37(2.09)
SAGPOOL | 68.12(3.07) | 69.96(2.48) | 67.51(2.65) | 49.91(2.22) | 51.07(2.31) | 67.23(2.15) | 68.83(1.13) | 67.51(2.51)
DIFFPOOL | 72.06(2.28) | 74.05(1.90) | 73.07(2.42) | 53.37(1.88) | 54.28(2.14) | 72.79(1.66) | 71.55(2.15) | 70.83(2.01)
BrainCheby | 75.08(1.98) | 76.14(2.38) | 74.09(1.84) | 56.25(2.12) | 57.37(2.05) | 72.55(2.45) | 73.36(1.88) | 72.62(1.33)
BrainNet-CNN | 74.09(2.49) | 73.71(1.96) | 73.27(2.21) | 54.03(2.20) | 55.25(2.46) | 68.37(1.71) | 69.97(1.30) | 68.51(2.02)
Ours w/o Contra. | 78.86(2.18) | 80.06(1.33) | 77.52(1.69) | 61.05(1.70) | 63.24(2.51) | 76.26(2.32) | 75.42(1.62) | 76.80(1.72)
Ours | 81.51(1.14) | 82.37(1.95) | 80.69(2.03) | 63.33(2.06) | 64.51(1.74) | 77.51(1.84) | 78.83(1.78) | 78.28(1.95)

####  V-D1 Experiment Setup

Six baseline models are utilized for comparison, including two machine
learning graph embedding models (t-BNE [56] and mCCA-ICA [57]), two deep graph
representation learning models designed for brain network embedding (BrainChey
[15] and BrainNet-CNN [14]), and two hierarchical graph neural networks with
graph pooling strategies (DIFFPOOL [20] and SAGPOOL [24]). Meanwhile, we
compare our model with and without optimizing contrastive loss to show that
the contrastive learning is beyond the data augmentation. The results for
gender and Alzheimer Disease (AD) classification are reported in accuracy,
precision and F1-score with their standard deviation (_std_). The results for
zygosity classification (i.e., $3$ classes classification task with class
labels as: not twins, monozygotic twins and dizygotic twins) are reported in
accuracy and Macro-F1-score with their _std_. The number of BUE and HGP
modules are set to $3$. We search the loss weights $\mu_{1}$ and $\mu_{2}$ in
range of $[0.1,1,5]$ and $[0.01,0.1,0.5,1]$ respectively and determine the
loss weights as $\mu_{1}=1$, $\mu_{2}=0.1$. The temperature parameter in
contrastive loss is set as $0.2$. Details of the hyperparameters analysis are
shown in section V-F.

####  V-D2 Results

Table I shows the results of gender classification, zygosity classification
and AD classification. It shows that our model achieves the best performance
comparing to all baseline methods on three tasks. For example, in the gender
classification, our model outperforms the baselines with at least $8.56\%$,
$8.18\%$ and $8.91\%$ increases in accuracy, precision and F1 scores,
respectively. In general, the deep graph neural networks are superior than the
traditional graph embedding methods (i.e., t-BNE and mCCA-ICA). When we remove
the supervision of the contrastive loss, the performance, though comparable to
baselines, decreases in a way. This manifests the effectiveness of the
contrastive learning and indicates that the contrastive learning is beyond a
data augmentation strategy.

Figure 3: The model performance obtained with different contrastive samples
generated by different window sizes. (A) shows the analysis on classification
tasks and (B) shows the analysis on regression tasks.

Figure 4: Loss weights analysis on classification tasks. (A) shows the
analysis on gender classification, (B) shows the analysis on zygosity
classification and (C) shows the analysis on AD classification. The red points
represent the best results. Figure 5: Brain saliency maps for classification
tasks. Here we identify: (1) top 15 regions associated with AD and NC from
OASIS, (2) top 10 regions associated with each sex and each zygosity from HCP.
Figure 6: Brain saliency maps for regression tasks. Here we identify: (1) top
15 regions associated with MMSE from OASIS, (2) top 10 regions associated with
Flanker score, Card-Sort score, Aggressive score, Intrusive score and Rule-
Break score from HCP.

###  V-E Regression tasks

Table II: Regression Mean Absolute Error (MAE) with s.t.d under 5-fold cross-validation. The values in bold show the best results. Method | OASIS | HCP
---|---|---
MMSE | Flanker | Card-Sort | Aggressive | Intrusive | Rule-Break
t-BNE | 2.02(0.36) | 1.69(0.19) | 1.58(0.22) | 1.89(0.10) | 1.84(0.22) | 1.77(0.41)
mCCA-ICA | 2.68(0.19) | 1.82(0.21) | 1.67(0.17) | 1.47(0.26) | 1.97(0.13) | 1.61(0.29)
SAGPOOL | 1.84(0.33) | 1.55(0.06) | 1.44(0.13) | 1.52(0.18) | 1.50(0.24) | 1.74(0.23)
DIFFPOOL | 1.27(0.20) | 1.34(0.14) | 1.16(0.30) | 1.27(0.41) | 1.25(0.07) | 1.43(0.15)
Brain-Cheby | 1.51(0.67) | 1.17(0.26) | 1.24(0.31) | 0.79(0.06) | 1.09(0.21) | 1.58(0.41)
BrainNetCNN | 1.26(0.19) | 1.43(0.24) | 0.91(0.11) | 1.33(0.23) | 1.14(0.13) | 1.29(0.19)
Ours w/o Contra. | 1.02(0.11) | 0.89(0.13) | 0.97(0.20) | 0.74(0.17) | 0.96(0.15) | 1.15(0.11)
Ours | 0.83(0.24) | 0.66(0.17) | 0.69(0.14) | 0.45(0.12) | 0.73(0.08) | 1.02(0.16)

####  V-E1 Experiment Setup

In the regression tasks, we use the same baselines for comparisons. The
regression tasks include predicting MMSE scores on OASIS data, Flanker scores,
Card-Sort scores, and $3$ ASR scores (i.e., Aggressive, Intrusive and Rule-
Break scores) on HCP data. Particularly, MMSE (Mini-Mental State Exam) test
[76], Flanker test [77] and Wisconsin Card-Sort test [78, 79, 80] are $3$
neuropsychological tests designed to measure the status and risks of human
neurodegenerative disease and mental illness. The ASR (Achenbach Adult Self-
Report) is a life function which is used to measure the emotion and social
support of adults. The structure of proposed model remains unchanged. The loss
weights are set as $\eta_{1}=0.5$ and $\eta_{2}=1$. The regression results are
reported in average Mean Absolute Errors (MAE) with its _std_ under 5-fold
cross validations.

####  V-E2 Results

The regression results are presented in Table II. It shows that our model
achieves the best MAE values comparing to all baseline methods. Similar to the
classification tasks, the deep graph neural networks are superior than
traditional graph embedding methods (i.e., t-BNE and mCCA-ICA). Comparing our
method with and without the supervision of the contrastive loss, we can hold
the conclusion that the contrastive learning can further boost the model
performance.

###  V-F Ablation Studies

In this section, we analyze the effect of $2$ hyperparameters on our model
performance. The first parameter is the window size ($d$) which we used to
clamp the fMRI bold signals when generating the contrastive functional brain
networks. Particularly, we set the window size as $[0,5,10,20,30,40,50]$,
respectively and generate different contrastive samples as the input of our
proposed model. Fig. 3 shows the performance of our framework under different
window sizes. It indicates that the best window size is around $d=10$. When
the window size decreases to $0$, the model performance declines since the
data is only duplicated without any substantial new samples. It is interesting
that the performance when $d=0$ is even worse than that obtained without
contrastive learning but with contrastive samples generated with $d=10$ (see
Ours w/o Contra. in Table I and II). The reason is that data augmentation is
introduced in the latter case, however, no augmented data is involved in the
first case.

We also analyze the effect of loss weights $\mu_{1}$ and $\mu_{2}$ on our
model performance. Fig. 4 presents the loss weight analysis on the $3$
classification tasks and the best results are achieved when $\mu_{1}=1$ and
$\mu_{2}=0.1$.

###  V-G Interpretation with Brain Saliency Map

We utilize the Class Activation Mapping (CAM) approach [81, 82, 83] to
generate the brain network saliency map, which indicates the top brain regions
associated with each prediction task. Figs 5 and 6 illustrate Brain Saliency
Maps for classification and regression tasks, respectively. For example, in
the classification task (AD vs. NC), the saliency map for AD highlights
multiple regions (such as Planum Polare, Frontal Operculum cortex,
Supracalcarine Cortex, etc.) which are conventionally conceived as the
biomarkers of AD in medical imaging analysis[84, 85, 86, 87]. In the meantime,
the saliency map for NC highlights many regions in Cerebellum and Frontal
lobe. These regions control cognitive thinking, motor control, and social
mentalizing as well as emotional self-experiences[88, 89, 90], in which AD
patients typically show problems. The details for all highlighted brain
regions for each task are summarized in the Supplementary Material. These
regions highlighted in the salience map can help us locating brain regions
associated with any phenotype, which deserve further clinical investigations.

##  VI Conclusion

We propose a novel contrastive learning framework with an interpretable
hierarchical signed graph representation learning model for brain functional
network mining. Additionally, a new data augmentation strategy is designed to
generate the contrastive samples for brain functional network data. Our new
framework is capable of generating more accurate representations for brain
functional networks in compared with other state-of-the-art methods and these
network representations can be used in various prediction tasks (e.g.,
classification and regression). Moreover, Brain saliency maps may assist with
phenotypic biomarker identification and provide interpretable explanation on
framework outcomes.

##  VII Acknowledgement

This study is partially supported by The National Institutes of Health
(R01AG071243, R01MH125928 and U01AG068057) and National Science Foundation
(IIS 2045848 and IIS 1837956).

Data were provided [in part] by the Human Connectome Project, MGH-USC
Consortium (Principal Investigators: Bruce R. Rosen, Arthur W. Toga and Van
Wedeen; U01MH093765) funded by the NIH Blueprint Initiative for Neuroscience
Research grant; the National Institutes of Health grant P41EB015896, and the
Instrumentation Grants S10RR023043, 1S10RR023401, 1S10RR019307.

Table III: The ROI names of the highlighted brain regions in the saliency map
in classification tasks

OASIS NC | OASIS AD | HCP Male | HCP Female | HCP Not Twins | HCP Monozygotic | HCP Dizygotic
---|---|---|---|---|---|---
| Paracingulate
---
Gyrus Right
Planum Polare Left | ctx-lh-precuneus | ctx-rh-superiorfrontal |  | ctx-lh-
---
lateraloccipital
| ctx-lh-
---
isthmuscingulate
ctx-lh-postcentral
| Paracingulate
---
Gyrus Right
| Frontal Operculum
---
Cortex Left
| ctx-rh-
---
superiorparietal
Right-Accumbens-area | ctx-rh-bankssts | ctx-rh-pericalcarine |  | ctx-lh-
---
caudalanteriorcingulate
Frontal Pole Right |  | Supracalcarine
---
Cortex Left
Right-Hippocampus |  | ctx-rh-
---
caudalmiddlefrontal
ctx-rh-parsorbitalis | ctx-rh-parsorbitalis | ctx-rh-parsorbitalis
Cerebellum 6 Right |  | Superior Temporal
---
Gyrus, anterior
division Right
| ctx-rh-
---
parahippocampal
ctx-lh-parsorbitalis | ctx-lh-precentral | ctx-rh-frontalpole | Right-Putamen
| Paracingulate
---
Gyrus Left
| Supramarginal
---
Gyrus, anterior
division Right
Right-Amygdala | Right-Amygdala |  | ctx-lh-
---
parahippocampal
ctx-lh-fusiform | ctx-rh-precentral
Left-Putamen | Left-Caudate | ctx-lh-pericalcarine | ctx-rh-paracentral | ctx-lh-entorhinal | ctx-lh-entorhinal |  | ctx-rh-
---
caudalmiddlefrontal
Cerebellum 8 Left |  | Middle Temporal
---
Gyrus, posterior
division Left
| ctx-lh-
---
transversetemporal
ctx-lh-precentral | ctx-lh-superiorfrontal |  | ctx-lh-
---
superiorfrontal
ctx-lh-precuneus
Cerebellum 7b Right |  | Superior Temporal
---
Gyrus, posterior
division Left
| ctx-rh-
---
transversetemporal
| ctx-lh-
---
isthmuscingulate
Right-Pallidum |  | ctx-lh-
---
temporalpole
ctx-lh-temporalpole
Heschl’s Gyrus Left | Heschl’s Gyrus Left |  | ctx-rh-
---
lateralorbitofrontal
| ctx-rh-
---
isthmuscingulate
| ctx-lh-
---
superiortemporal
| ctx-lh-
---
superiorparietal
| ctx-rh-
---
transversetemporal
Cuneal Cortex Right |  | Intracalcarine
---
Cortex Left
| ctx-lh-
---
temporalpole
| ctx-lh-
---
caudalanteriorcingulate
| ctx-rh-
---
caudalmiddlefrontal
Left-Pallidum |  | ctx-rh-
---
transversetemporal
| Lateral Occipital
---
Cortex,
superior division Right
| Middle Frontal
---
Gyrus Left
|  |  |  |
Precuneous Cortex | Planum Polare Right |  |  |  |  |
Cerebellum Crus2 Left |  | Temporal Fusiform
---
Cortex, anterior
division Left
|  |  |  |
Brain-Stem |  | Middle Temporal
---
Gyrus, temporooccipital
part Left
|  |  |  |
Cerebellum 8 Right |  | Supracalcarine
---
Cortex Right
|  |  |  |

Table IV: The ROI names of the highlighted brain regions in the saliency map in regression tasks OASIS MMSE | Flanker | Card-Sort | Aggressive | Intrusive | Rule-Break
---|---|---|---|---|---
Right-Caudate | Left-Accumbens-area | Left-Accumbens-area | ctx-lh-bankssts | ctx-lh-bankssts | ctx-lh-precuneus
Temporal Pole Right | ctx-rh-fusiform | Left-Putamen |  | ctx-lh-
---
middletemporal
ctx-lh-parsorbitalis | ctx-lh-lingual
| Middle Temporal
---
Gyrus, posterior
division Right
ctx-lh-inferiortemporal |  | ctx-lh-
---
caudalmiddlefrontal
| ctx-lh-
---
inferiortemporal
| ctx-lh-
---
inferiortemporal
| ctx-lh-
---
inferiortemporal
Cerebellum Crus1 Right | ctx-rh-insula | ctx-rh-frontalpole |  | ctx-lh-
---
lateraloccipital
| ctx-lh-
---
parahippocampal
Right-Caudate
| Temporal Occipital
---
Fusiform Cortex Left
| ctx-lh-
---
middletemporal
| ctx-lh-
---
rostralanteriorcingulate
ctx-lh-precentral |  | ctx-lh-
---
caudalanteriorcingulate
| ctx-rh-
---
lateraloccipital
Planum Temporale Left | ctx-lh-postcentral |  | ctx-rh-
---
caudalmiddlefrontal
ctx-rh-temporalpole |  | ctx-rh-
---
caudalmiddlefrontal
| ctx-rh-
---
temporalpole
| Middle Temporal
---
Gyrus, temporooccipital
part Left
ctx-lh-frontalpole |  | ctx-rh-
---
middletemporal
ctx-rh-frontalpole |  | ctx-rh-
---
supramarginal
| ctx-lh-
---
supramarginal
| Temporal Occipital
---
Fusiform Cortex Right
ctx-lh-temporalpole | ctx-lh-frontalpole | ctx-rh-parsorbitalis | ctx-rh-paracentral | ctx-rh-insula
Planum Temporale Right |  | ctx-rh-
---
superiortemporal
ctx-rh-precentral |  | ctx-rh-
---
parstriangularis
| ctx-rh-
---
parstriangularis
| ctx-rh-
---
parstriangularis
Frontal Orbital Cortex Left | ctx-rh-precentral | ctx-rh-precuneus | ctx-rh-entorhinal | ctx-lh-precentral | Right-Amygdala
| Middle Temporal
---
Gyrus, posterior
division Left
|  |  |  |
Vermis 9 |  |  |  |  |
| Middle Temporal
---
Gyrus, temporooccipital
part Right
|  |  |  |
Left-Caudate |  |  |  |  |
Temporal Pole Left |  |  |  |  |

## References

  * [1] H. Rusinek, S. De Santi, D. Frid, W.-H. Tsui, C. Y. Tarshish, A. Convit, and M. J. de Leon, “Regional brain atrophy rate predicts future cognitive decline: 6-year longitudinal mr imaging study of normal aging,” _Radiology_ , vol. 229, no. 3, pp. 691–696, 2003.
  * [2] M. R. Sabuncu and E. Konukoglu, “Clinical prediction from structural brain mri scans: a large-scale empirical study,” _Neuroinformatics_ , vol. 13, no. 1, pp. 31–46, 2015.
  * [3] S. Seo, J. Mohr, A. Beck, T. Wüstenberg, A. Heinz, and K. Obermayer, “Predicting the future relapse of alcohol-dependent patients from structural and functional brain images,” _Addiction biology_ , vol. 20, no. 6, pp. 1042–1055, 2015.
  * [4] M. P. Van Den Heuvel, R. S. Kahn, J. Goñi, and O. Sporns, “High-cost, high-capacity backbone for global brain communication,” _Proceedings of the National Academy of Sciences_ , vol. 109, no. 28, pp. 11 372–11 377, 2012\.
  * [5] O. Sporns, “The human connectome: origins and challenges,” _Neuroimage_ , vol. 80, pp. 53–61, 2013.
  * [6] M. G. Mattar and D. S. Bassett, “Brain network architecture,” _Network science in cognitive psychology_ , p. 30, 2019.
  * [7] Y. Zhang, L. Zhan, S. Wu, P. Thompson, and H. Huang, “Disentangled and proportional representation learning for multi-view brain connectomes,” in _International Conference on Medical Image Computing and Computer-Assisted Intervention_. Springer, 2021, pp. 508–518.
  * [8] R. E. Beaty, Y. N. Kenett, A. P. Christensen, M. D. Rosenberg, M. Benedek, Q. Chen, A. Fink, J. Qiu, T. R. Kwapil, M. J. Kane _et al._ , “Robust prediction of individual creative ability from brain functional connectivity,” _Proceedings of the National Academy of Sciences_ , vol. 115, no. 5, pp. 1087–1092, 2018.
  * [9] C. J. Brown, K. P. Moriarty, S. P. Miller, B. G. Booth, J. G. Zwicker, R. E. Grunau, A. R. Synnes, V. Chau, and G. Hamarneh, “Prediction of brain network age and factors of delayed maturation in very preterm infants,” in _International Conference on Medical Image Computing and Computer-Assisted Intervention_. Springer, 2017, pp. 84–91.
  * [10] T. Eichele, S. Debener, V. D. Calhoun, K. Specht, A. K. Engel, K. Hugdahl, D. Y. Von Cramon, and M. Ullsperger, “Prediction of human errors by maladaptive changes in event-related brain networks,” _Proceedings of the National Academy of Sciences_ , vol. 105, no. 16, pp. 6173–6178, 2008.
  * [11] X. Li, Y. Li, and X. Li, “Predicting clinical outcomes of alzheimer’s disease from complex brain networks,” in _International Conference on Advanced Data Mining and Applications_. Springer, 2017, pp. 519–525.
  * [12] D. E. Warren, N. L. Denburg, J. D. Power, J. Bruss, E. J. Waldron, H. Sun, S. E. Petersen, and D. Tranel, “Brain network theory can predict whether neuropsychological outcomes will differ from clinical expectations,” _Archives of Clinical Neuropsychology_ , vol. 32, no. 1, pp. 40–52, 2017\.
  * [13] C. Hu, R. Ju, Y. Shen, P. Zhou, and Q. Li, “Clinical decision support for alzheimer’s disease based on deep learning and brain network,” in _2016 IEEE International Conference on Communications (ICC)_. IEEE, 2016, pp. 1–6.
  * [14] J. Kawahara, C. J. Brown, S. P. Miller, B. G. Booth, V. Chau, R. E. Grunau, J. G. Zwicker, and G. Hamarneh, “Brainnetcnn: Convolutional neural networks for brain networks; towards predicting neurodevelopment,” _NeuroImage_ , vol. 146, pp. 1038–1049, 2017.
  * [15] S. I. Ktena, S. Parisot, E. Ferrante, M. Rajchl, M. Lee, B. Glocker, and D. Rueckert, “Metric learning with spectral graph convolutions on brain connectivity networks,” _NeuroImage_ , vol. 169, pp. 431–442, 2018.
  * [16] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” _arXiv preprint arXiv:1609.02907_ , 2016.
  * [17] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Lio, and Y. Bengio, “Graph attention networks,” _arXiv preprint arXiv:1710.10903_ , 2017.
  * [18] Z. Ying, D. Bourgeois, J. You, M. Zitnik, and J. Leskovec, “Gnnexplainer: Generating explanations for graph neural networks,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [19] Y. Zhang, L. Zhan, W. Cai, P. Thompson, and H. Huang, “Integrating heterogeneous brain networks for predicting brain disease conditions,” in _International Conference on Medical Image Computing and Computer-Assisted Intervention_. Springer, 2019, pp. 214–222.
  * [20] Z. Ying, J. You, C. Morris, X. Ren, W. Hamilton, and J. Leskovec, “Hierarchical graph representation learning with differentiable pooling,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [21] C. C. Hilgetag and A. Goulas, “‘hierarchy’in the organization of brain networks,” _Philosophical Transactions of the Royal Society B_ , vol. 375, no. 1796, p. 20190319, 2020.
  * [22] R. Mastrandrea, A. Gabrielli, F. Piras, G. Spalletta, G. Caldarelli, and T. Gili, “Organization and hierarchy of the human functional brain network lead to a chain-like core,” _Scientific reports_ , vol. 7, no. 1, pp. 1–13, 2017.
  * [23] D. Meunier, R. Lambiotte, A. Fornito, K. Ersche, and E. T. Bullmore, “Hierarchical modularity in human brain functional networks,” _Frontiers in neuroinformatics_ , vol. 3, p. 37, 2009.
  * [24] J. Lee, I. Lee, and J. Kang, “Self-attention graph pooling,” in _International conference on machine learning_. PMLR, 2019, pp. 3734–3743.
  * [25] Z. Zhang, J. Bu, M. Ester, J. Zhang, C. Yao, Z. Yu, and C. Wang, “Hierarchical graph pooling with structure learning,” _arXiv preprint arXiv:1911.05954_ , 2019.
  * [26] X. Li, Y. Zhou, N. Dvornek, M. Zhang, S. Gao, J. Zhuang, D. Scheinost, L. H. Staib, P. Ventola, and J. S. Duncan, “Braingnn: Interpretable brain graph neural network for fmri analysis,” _Medical Image Analysis_ , vol. 74, p. 102233, 2021.
  * [27] H. Tang, G. Ma, L. He, H. Huang, and L. Zhan, “Commpool: An interpretable graph pooling framework for hierarchical graph representation learning,” _Neural Networks_ , vol. 143, pp. 669–677, 2021.
  * [28] D. Cartwright and F. Harary, “Structural balance: a generalization of heider’s theory.” _Psychological review_ , vol. 63, no. 5, p. 277, 1956.
  * [29] T. Derr, Y. Ma, and J. Tang, “Signed graph convolutional networks,” in _2018 IEEE International Conference on Data Mining (ICDM)_. IEEE, 2018, pp. 929–934.
  * [30] F. Heider, “Attitudes and cognitive organization,” _The Journal of psychology_ , vol. 21, no. 1, pp. 107–112, 1946.
  * [31] Y. Li, Y. Tian, J. Zhang, and Y. Chang, “Learning signed network embedding via graph attention,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 04, 2020, pp. 4772–4779.
  * [32] P. Khosla, P. Teterwak, C. Wang, A. Sarna, Y. Tian, P. Isola, A. Maschinot, C. Liu, and D. Krishnan, “Supervised contrastive learning,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 18 661–18 673, 2020\.
  * [33] D. Berthelot, N. Carlini, I. Goodfellow, N. Papernot, A. Oliver, and C. A. Raffel, “Mixmatch: A holistic approach to semi-supervised learning,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [34] Q. Xie, Z. Dai, E. Hovy, T. Luong, and Q. Le, “Unsupervised data augmentation for consistency training,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 6256–6268, 2020.
  * [35] K. Hassani and A. H. Khasahmadi, “Contrastive multi-view representation learning on graphs,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 4116–4126.
  * [36] Y. You, T. Chen, Y. Sui, T. Chen, Z. Wang, and Y. Shen, “Graph contrastive learning with augmentations,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 5812–5823, 2020.
  * [37] Y. Zhu, Y. Xu, F. Yu, Q. Liu, S. Wu, and L. Wang, “Graph contrastive learning with adaptive augmentation,” in _Proceedings of the Web Conference 2021_ , 2021, pp. 2069–2080.
  * [38] T. Zhao, Y. Liu, L. Neves, O. Woodford, M. Jiang, and N. Shah, “Data augmentation for graph neural networks,” _arXiv preprint arXiv:2006.06830_ , 2020.
  * [39] M. P. van den Heuvel and O. Sporns, “Network hubs in the human brain,” _Trends in cognitive sciences_ , vol. 17, no. 12, pp. 683–696, 2013.
  * [40] M. U. Ilyas, M. Z. Shafiq, A. X. Liu, and H. Radha, “A distributed and privacy preserving algorithm for identifying information hubs in social networks,” in _2011 Proceedings IEEE INFOCOM_. IEEE, 2011, pp. 561–565.
  * [41] K. Hwang, M. N. Hallquist, and B. Luna, “The development of hub architecture in the human functional brain network,” _Cerebral Cortex_ , vol. 23, no. 10, pp. 2380–2393, 2013.
  * [42] J. Chen, T. Ma, and C. Xiao, “Fastgcn: fast learning with graph convolutional networks via importance sampling,” _arXiv preprint arXiv:1801.10247_ , 2018\.
  * [43] W. Huang, T. Zhang, Y. Rong, and J. Huang, “Adaptive sampling towards fast graph representation learning,” in _Advances in neural information processing systems_ , 2018, pp. 4558–4567.
  * [44] H. Dai, B. Dai, and L. Song, “Discriminative embeddings of latent variable models for structured data,” in _International conference on machine learning_ , 2016, pp. 2702–2711.
  * [45] D. K. Duvenaud, D. Maclaurin, J. Iparraguirre, R. Bombarell, T. Hirzel, A. Aspuru-Guzik, and R. P. Adams, “Convolutional networks on graphs for learning molecular fingerprints,” in _Advances in neural information processing systems_ , 2015, pp. 2224–2232.
  * [46] J. Liu, G. Ma, F. Jiang, C.-T. Lu, S. Y. Philip, and A. B. Ragin, “Community-preserving graph convolutions for structural and functional joint embedding of brain networks,” in _2019 IEEE International Conference on Big Data (Big Data)_. IEEE, 2019, pp. 1163–1168.
  * [47] W. Hamilton, Z. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” in _Advances in neural information processing systems_ , 2017, pp. 1024–1034.
  * [48] Y. Li, D. Tarlow, M. Brockschmidt, and R. Zemel, “Gated graph sequence neural networks,” _arXiv preprint arXiv:1511.05493_ , 2015.
  * [49] O. Vinyals, S. Bengio, and M. Kudlur, “Order matters: Sequence to sequence for sets,” _arXiv preprint arXiv:1511.06391_ , 2015.
  * [50] M. Zhang, Z. Cui, M. Neumann, and Y. Chen, “An end-to-end deep learning architecture for graph classification,” in _Thirty-Second AAAI Conference on Artificial Intelligence_ , 2018.
  * [51] G. Ma, N. K. Ahmed, T. L. Willke, and P. S. Yu, “Deep graph similarity learning: A survey,” _Data Mining and Knowledge Discovery_ , vol. 35, no. 3, pp. 688–725, 2021.
  * [52] H. Gao and S. Ji, “Graph u-nets,” in _international conference on machine learning_. PMLR, 2019, pp. 2083–2092.
  * [53] H. Yuan and S. Ji, “Structpool: Structured graph pooling via conditional random fields,” in _Proceedings of the 8th International Conference on Learning Representations_ , 2020.
  * [54] L. E. Korthauer, L. Zhan, O. Ajilore, A. Leow, and I. Driscoll, “Disrupted topology of the resting state structural connectome in middle-aged apoe $\varepsilon$4 carriers,” _Neuroimage_ , vol. 178, pp. 295–305, 2018.
  * [55] L. Zhan, Y. Liu, J. Zhou, J. Ye, and P. M. Thompson, “Boosting classification accuracy of diffusion mri derived brain networks for the subtypes of mild cognitive impairment using higher order singular value decomposition,” in _2015 IEEE 12th International Symposium on Biomedical Imaging (ISBI)_. IEEE, 2015, pp. 131–135.
  * [56] B. Cao, L. He, X. Wei, M. Xing, P. S. Yu, H. Klumpp, and A. D. Leow, “t-bne: Tensor-based brain network embedding,” in _Proceedings of the 2017 SIAM International Conference on Data Mining_. SIAM, 2017, pp. 189–197.
  * [57] J. Sui, G. Pearlson, A. Caprihan, T. Adali, K. A. Kiehl, J. Liu, J. Yamamoto, and V. D. Calhoun, “Discriminating schizophrenia and bipolar disorder by fusing fmri and dti in a multimodal cca+ joint ica model,” _Neuroimage_ , vol. 57, no. 3, pp. 839–855, 2011.
  * [58] Y. Zhang and H. Huang, “New graph-blind convolutional network for brain connectome data analysis,” in _International Conference on Information Processing in Medical Imaging_. Springer, 2019, pp. 669–681.
  * [59] H. Jiang, P. Cao, M. Xu, J. Yang, and O. Zaiane, “Hi-gcn: A hierarchical graph convolution network for graph embedding learning of brain network and brain disorders prediction,” _Computers in Biology and Medicine_ , vol. 127, p. 104096, 2020.
  * [60] J. Jung, J. Yoo, and U. Kang, “Signed graph diffusion network,” _arXiv preprint arXiv:2012.14191_ , 2020.
  * [61] X. Shen and F.-L. Chung, “Deep network embedding for graph representation learning in signed networks,” _IEEE transactions on cybernetics_ , vol. 50, no. 4, pp. 1556–1568, 2018.
  * [62] Q. Huang, M. Yamada, Y. Tian, D. Singh, D. Yin, and Y. Chang, “Graphlime: Local interpretable model explanations for graph neural networks,” _arXiv preprint arXiv:2001.06216_ , 2020.
  * [63] H. Cui, W. Dai, Y. Zhu, X. Li, L. He, and C. Yang, “Brainnnexplainer: An interpretable graph neural network framework for brain network based disease analysis,” _arXiv preprint arXiv:2107.05097_ , 2021.
  * [64] D. Xu, W. Cheng, D. Luo, H. Chen, and X. Zhang, “Infogcl: Information-aware graph contrastive learning,” _Advances in Neural Information Processing Systems_ , vol. 34, 2021.
  * [65] Z. Wu, Y. Xiong, S. X. Yu, and D. Lin, “Unsupervised feature learning via non-parametric instance discrimination,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 3733–3742.
  * [66] K. Sohn, “Improved deep metric learning with multi-class n-pair loss objective,” _Advances in neural information processing systems_ , vol. 29, 2016.
  * [67] A. Van den Oord, Y. Li, and O. Vinyals, “Representation learning with contrastive predictive coding,” _arXiv e-prints_ , pp. arXiv–1807, 2018\.
  * [68] D. C. Van Essen, S. M. Smith, D. M. Barch, T. E. Behrens, E. Yacoub, K. Ugurbil, W.-M. H. Consortium _et al._ , “The wu-minn human connectome project: an overview,” _Neuroimage_ , vol. 80, pp. 62–79, 2013.
  * [69] P. J. LaMontagne, T. L. Benzinger, J. C. Morris, S. Keefe, R. Hornbeck, C. Xiong, E. Grant, J. Hassenstab, K. Moulder, A. G. Vlassenko _et al._ , “Oasis-3: longitudinal neuroimaging, clinical, and cognitive dataset for normal aging and alzheimer disease,” _MedRxiv_ , 2019.
  * [70] S. Whitfield-Gabrieli and A. Nieto-Castanon, “Conn: a functional connectivity toolbox for correlated and anticorrelated brain networks,” _Brain connectivity_ , vol. 2, no. 3, pp. 125–141, 2012.
  * [71] I. Fortel, L. E. Korthauer, Z. Morrissey, L. Zhan, O. Ajilore, O. Wolfson, I. Driscoll, D. Schonfeld, and A. Leow, “Connectome signatures of hyperexcitation in cognitively intact middle-aged female apoe-$\varepsilon$4 carriers,” _Cerebral Cortex_ , vol. 30, no. 12, pp. 6350–6362, 2020.
  * [72] O. Ajilore, L. Zhan, J. GadElkarim, A. Zhang, J. Feusner, S. Yang, P. M. Thompson, A. Kumar, and A. Leow, “Constructing the resting state structural connectome,” _Frontiers in neuroinformatics_ , vol. 7, p. 30, 2013.
  * [73] B. Fischl, “Freesurfer,” _Neuroimage_ , vol. 62, no. 2, pp. 774–781, 2012\.
  * [74] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [75] O. Shchur, M. Mumme, A. Bojchevski, and S. Günnemann, “Pitfalls of graph neural network evaluation,” _arXiv preprint arXiv:1811.05868_ , 2018.
  * [76] T. N. Tombaugh and N. J. McIntyre, “The mini-mental state examination: a comprehensive review,” _Journal of the American Geriatrics Society_ , vol. 40, no. 9, pp. 922–935, 1992.
  * [77] B. A. Eriksen and C. W. Eriksen, “Effects of noise letters upon the identification of a target letter in a nonsearch task,” _Perception & psychophysics_, vol. 16, no. 1, pp. 143–149, 1974.
  * [78] V. C. Pangman, J. Sloan, and L. Guse, “An examination of psychometric properties of the mini-mental state examination and the standardized mini-mental state examination: implications for clinical practice,” _Applied Nursing Research_ , vol. 13, no. 4, pp. 209–213, 2000.
  * [79] O. Monchi, M. Petrides, V. Petre, K. Worsley, and A. Dagher, “Wisconsin card sorting revisited: distinct neural circuits participating in different stages of the task identified by event-related functional magnetic resonance imaging,” _Journal of Neuroscience_ , vol. 21, no. 19, pp. 7733–7741, 2001\.
  * [80] E. A. Berg, “A simple objective technique for measuring flexibility in thinking,” _The Journal of general psychology_ , vol. 39, no. 1, pp. 15–22, 1948.
  * [81] W. Zhang, L. Zhan, P. Thompson, and Y. Wang, “Deep representation learning for multimodal brain networks,” in _International Conference on Medical Image Computing and Computer-Assisted Intervention_. Springer, 2020, pp. 613–624.
  * [82] S. Arslan, S. I. Ktena, B. Glocker, and D. Rueckert, “Graph saliency maps through spectral convolutional networks: Application to sex classification with brain connectivity,” in _Graphs in biomedical image analysis and integrating medical imaging and non-imaging modalities_. Springer, 2018, pp. 3–13.
  * [83] P. E. Pope, S. Kolouri, M. Rostami, C. E. Martin, and H. Hoffmann, “Explainability methods for graph convolutional neural networks,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 10 772–10 781.
  * [84] J. Rasero, N. Amoroso, M. La Rocca, S. Tangaro, R. Bellotti, S. Stramaglia, and A. D. N. Initiative, “Multivariate regression analysis of structural mri connectivity matrices in alzheimer’s disease,” _PLoS One_ , vol. 12, no. 11, p. e0187281, 2017.
  * [85] M. Kutova, J. Mrzilkova, J. Riedlova, and P. Zach, “Asymmetric changes in limbic cortex and planum temporale in patients with alzheimer disease,” _Current Alzheimer Research_ , vol. 15, no. 14, pp. 1361–1368, 2018.
  * [86] L. V. Hiscox, C. L. Johnson, M. D. McGarry, H. Marshall, C. W. Ritchie, E. J. Van Beek, N. Roberts, and J. M. Starr, “Mechanical property alterations across the cerebral cortex due to alzheimer’s disease,” _Brain communications_ , vol. 2, no. 1, p. fcz049, 2020.
  * [87] A. Hafkemeijer, C. Möller, E. G. Dopper, L. C. Jiskoot, T. M. Schouten, J. C. Van Swieten, W. M. Van der Flier, H. Vrenken, Y. A. Pijnenburg, F. Barkhof _et al._ , “Resting state functional connectivity differences between behavioral variant frontotemporal dementia and alzheimer’s disease,” _Frontiers in human neuroscience_ , vol. 9, p. 474, 2015.
  * [88] C. J. Stoodley, E. M. Valera, and J. D. Schmahmann, “Functional topography of the cerebellum for motor and cognitive tasks: an fmri study,” _Neuroimage_ , vol. 59, no. 2, pp. 1560–1570, 2012.
  * [89] F. Van Overwalle, Q. Ma, and E. Heleven, “The posterior crus ii cerebellum is specialized for social mentalizing and emotional self-experiences: a meta-analysis,” _Social Cognitive and Affective Neuroscience_ , vol. 15, no. 9, pp. 905–928, 2020.
  * [90] R. P. Sawyer, F. Rodriguez-Porcel, M. Hagen, R. Shatz, and A. J. Espay, “Diagnosing the frontal variant of alzheimer’s disease: a clinician’s yellow brick road,” _Journal of clinical movement disorders_ , vol. 4, no. 1, pp. 1–9, 2017.


# ESDIRK-based nonlinear model predictive control
for stochastic differential-algebraic equations

Anders Hilmar Damm Christensen, Nicola Cantisani, John Bagterp Jørgensen A. H.
D. Christensen, N. Cantisani, and J. B. Jørgensen are with the Department of
Applied Mathematics and Computer Science, Technical University of Denmark,
DK-2800 Kgs. Lyngby, Denmark. Corresponding author: J. B. Jørgensen (E-mail:
jbjo@dtu.dk<EMAIL_ADDRESS>
###### Abstract

In this paper, we present a nonlinear model predictive control (NMPC)
algorithm for systems modeled by semi-explicit stochastic differential-
algebraic equations (DAEs) of index 1. The NMPC combines a continuous-discrete
extended Kalman filter (CD-EKF) with an optimal control problem (OCP) for
setpoint tracking. We discretize the OCP using direct multiple shooting. We
apply an explicit singly diagonal implicit Runge-Kutta (ESDIRK) integration
scheme to solve systems of DAEs, both for the one-step prediction in the CD-
EKF and in each shooting interval of the discretized OCP. The ESDIRK methods
use an iterated internal numerical differentiation approach for precise
sensitivity computations. These sensitivities are used to provide accurate
gradient information in the OCP and to efficiently integrate the covariance
information in the CD-EKF. Subsequently, we present a simulation case study
where we apply the NMPC to a simple alkaline electrolyzer stack model. We use
the NMPC to track a time-varying setpoint for the stack temperature subject to
input bound constraints.

##  I Introduction

Nonlinear model predictive control (NMPC) algorithms are widely used in the
control of chemical systems to optimize economics or minimize deviations from
setpoints [1]. An NMPC algorithm is a feedback control strategy that applies
the solutions from optimal control problems (OCPs) based on a mathematical
model of the system in a moving horizon fashion. Since chemical systems are
best modeled using differential-algebraic equations (DAEs) [2], it follows
that the models employed in NMPC strategies should also be formulated as DAEs.

Implementing the OCPs using direct multiple shooting discretization requires
numerical integration schemes to solve initial-value problems (IVPs). Software
solutions for solving IVPs for DAEs include the codes DASSL, which utilizes a
backward differentiation formula (BDF), and RADAU5, employing a 3-stages fully
implicit Runge-Kutta (IRK) integration scheme [3]. Due to the frequent
discontinuities inherent in multiple shooting-based NMPC, BDF methods may
experience frequent restarts, leading to a reduction in integration order [4].
In contrast, the fully implicit IRK methods do not suffer from the reduction
of integration order. However, these methods are in general considered
computationally expensive methods as they require the solution of large
systems of nonlinear equations in each integration step. Therefore, the
diagonal IRK methods, such as explicit singly diagonal implicit Runge-Kutta
(ESDIRK) integration schemes, offer a computationally efficient alternative.
ESDIRK methods are characterized by their identical coefficients in the
diagonal of their Butcher tableaus with the first step being explicit. Such
methods enable the reuse of matrix factorizations across stages, thereby
enhancing computational efficiency in solving IVPs for DAEs.

Solving the OCPs also requires integrator sensitivity information. Equipping
the implementation of the ESDIRK integration schemes with iterated sensitivity
computations based on the principle of internal numerical differentiation
(IND), ensures precise sensitivity information. The principle of IND involves
computing the sensitivities by directly differentiating the discretization
scheme generated adaptively from the integrator [5, 6]. [7] presents iterated
IND sensitivity calculations for ESDIRK-based OCPs applied to ODEs. It is
demonstrated that such sensitivities can enhance convergence properties in an
NLP solver compared to approximate sensitivities.
Extended Kalman filters are widely used state estimators in NMPC applications
[8]. [9] presents a continuous-discrete extended Kalman filter (CD-EKF) for
systems described by stochastic DAEs (SDAEs). This CD-EKF implementation is
computationally efficient due to the use of ESDIRK methods for simultaneous
integration of the mean-covariance pair. The implementation obtains the
covariance information by using sensitivity computations of the ESDIRK method.
Here, the sensitivities are computed using staggered direct approach based on
the assumption of a constant Jacobian in each accepted integration step.

This paper presents an ESDIRK-based NMPC algorithm for setpoint tracking in
continuous-discrete semi-explicit index-1 SDAEs. The NMPC algorithm combines a
CD-EKF with an OCP discretized using direct multiple shooting. The novelty in
this paper is the direct demonstration of how the ESDIRK methods with iterated
IND sensitivity computations can be applied to address both the state
estimation problem and the OCP, resulting in a computationally efficient NMPC
solution for SDAEs. We test the NMPC algorithm on a numerical case study of an
alkaline electrolyzer stack affected by a stochastic process. We use the NMPC
to track a time-varying setpoint for the temperature in the stack.

The rest of the paper is structured as follows. Section II introduces the
SDAEs. In Section III, we describe the ESDIRK integration scheme for solving
the IVPs of DAEs in the CD-EKF and the OCP. We derive iterated IND
sensitivities for these ESDIRK methods. We also introduce an implicit-explicit
scheme for solving the IVPs for SDAE systems. In Section IV, we introduce the
NMPC and demonstrate how to implement the CD-EKF and OCP for index-1 semi-
explicit DAEs. Section V presents a case study using the model of an alkaline
electrolyzer stack. Section VI concludes the paper.

##  II Stochastic differential-algebraic equations

We consider the continuous-discrete semi-explicit SDAEs of the form

$\displaystyle dx(t)$ $\displaystyle=f(t,x(t),y(t),u(t),d(t))dt+\sigma
d\omega(t),$ (1a) $\displaystyle 0$ $\displaystyle=g(t,x(t),y(t),u(t),d(t)),$
(1b) $\displaystyle y_{k}^{m}$
$\displaystyle=m(t_{k},x_{k},y_{k},u_{k},d_{k})+v_{k},$ (1c) $\displaystyle
z(t)$ $\displaystyle=h(t,x(t),y(t),u(t),d(t)).$ (1d)

$t$ is time, $x(t)$ is the vector of differential states, $y(t)$ is the vector
of algebraic states, $u(t)$ is the vector of inputs, and $d(t)$ is the vector
of disturbances. $\omega(t)$ is a standard Wiener process, i.e.,
$d\omega(t)\sim N_{iid}(0,Idt)$, representing the plant-model mismatch and
$v_{k}\sim N_{iid}(0,R_{k})$ is the measurement noise with $R_{k}$ being the
covariance. $y_{k}^{m}=y^{m}(t_{k})$ represent the measurements at the time
$t_{k}$ and $z(t)$ is the output of the model. We assume the system in (1) is
of index-1, i.e., $(\partial g/\partial y)^{-1}$ exists.

##  III Integration and sensitivity generation

In this section, we introduce the ESDIRK integration schemes to solve IVPs for
DAEs of index-1. For these ESDIRK methods, we demonstrate how to generate the
forward sensitivities using the principle of IND.

TABLE I: Butcher tableaus for some ESDIRK methods.

ESDIRK12
---
$0$ | $0$ |
$1$ | $b_{1}$ | $\gamma$
$x_{k+1}$ | $b_{1}$ | $\gamma$
$\hat{x}_{k+1}$ | $\hat{b}_{1}$ | $\hat{b}_{2}$

ESDIRK23
---
$0$ | $0$ |  |
$c_{2}$ | $a_{21}$ | $\gamma$ |
$1$ | $b_{1}$ | $b_{1}$ | $\gamma$
$x_{k+1}$ | $b_{1}$ | $b_{1}$ | $\gamma$
$\hat{x}_{k+1}$ | $\hat{b}_{1}$ | $\hat{b}_{2}$ | $\hat{b}_{3}$

ESDIRK34
---
$0$ | $0$ |  |  |
$c_{2}$ | $a_{21}$ | $\gamma$ |  |
$c_{3}$ | $a_{31}$ | $a_{32}$ | $\gamma$ |
$1$ | $b_{1}$ | $b_{2}$ | $b_{3}$ | $\gamma$
$x_{k+1}$ | $b_{1}$ | $b_{2}$ | $b_{3}$ | $\gamma$
$\hat{x}_{k+1}$ | $\hat{b}_{1}$ | $\hat{b}_{2}$ | $\hat{b}_{3}$ | $\hat{b}_{4}$

###  III-A The ESDIRK integration schemes for DAEs

We solve the IVPs

$\displaystyle\frac{dx(t)}{dt}$ $\displaystyle=f(t,x(t),y(t),u,d),\quad
x(t_{0})=x_{0},$ (2a) $\displaystyle 0$
$\displaystyle=g(t,x(t),y(t),u,d),\quad y(t_{0})=y_{0},$ (2b)

for $t\in[t_{0},\,t_{f}]$ using the ESDIRK integration schemes. The
$n_{s}$-stages ESDIRK methods for solving (2a) are [10]

$\displaystyle T_{i}$ $\displaystyle=t_{k}+c_{i}h_{k},$ (3a) $\displaystyle
X_{i}$
$\displaystyle=x_{k}+h_{k}\sum_{j=1}^{i}a_{ij}f(T_{j},X_{j},Y_{j},u,d),$ (3b)
$\displaystyle x_{k+1}$
$\displaystyle=x_{k}+h_{k}\sum_{i=1}^{n_{s}}b_{i}f(T_{j},X_{j},Y_{j},u,d),$
(3c) $\displaystyle\hat{x}_{k+1}$
$\displaystyle=x_{k}+h_{k}\sum_{i=1}^{n_{s}}\hat{b}_{i}f(T_{j},X_{j},Y_{j},u,d).$
(3d)

$T_{i}$ are the internal nodes, $X_{i}$ are the differential stages, and
$Y_{i}$ algebraic stages at iteration $k$ for $i=1,\dots,n_{s}$. $x_{k}$ and
$x_{k+1}$ are the steps computed at $t_{k}$ and $t_{k+1}=t_{k}+h_{k}$. $h_{k}$
is the integration step size at iteration $k$ which may be controlled based on
a local error estimate, $e_{k+1}=x_{k+1}-\hat{x}_{k+1}$, using the embedded
method. To satisfy the consistency equations in (2b) we compute $Y_{j}$ such
that

$0=g(T_{j},X_{j},Y_{j},u,d),\quad\text{for}\quad j=1,\dots,n_{s}.$ (4)

Table I presents the Butcher tableaus for the ESDIRK12, ESDIRK23, and ESDIRK34
methods. The numerical values are provided by [11]. For the rest of this
section, we leave out the arguments $T_{j}$ and $d$.

We combine (3b) and (4) into the system of nonlinear equations

$R_{i}(S_{i})=\begin{bmatrix}X_{i}-h_{k}\gamma f(X_{i},Y_{i},u)-\psi_{i}\\\
-g(X_{i},Y_{i},u)\end{bmatrix}=0,$ (5)

for $i=2,\dots,n_{s}$. $S_{i}=\begin{bmatrix}X_{i};&Y_{i}\end{bmatrix}$ and

$\psi_{i}=x_{k}+h_{k}\sum_{j=1}^{i-1}a_{ij}f(X_{j},Y_{j},u),$ (6)

and solve (5) in each integration step. We do this using the inexact Newton
scheme

$\displaystyle M_{k}\Delta S_{i}^{[l]}$ $\displaystyle=-R_{i}(S_{i}^{[l]}),$
(7a) $\displaystyle S_{i}^{[l+1]}$ $\displaystyle=S_{i}^{[l]}+\Delta
S_{i}^{[l]},$ (7b)

for $l=0,\dots,v_{i,k}-1$. $v_{i,k}$ is the number of Newton-type iterations
required for stage $i$ at iteration $k$ to satisfy some convergence criteria.
We choose the convergence criteria

$\left\|R_{i}(S_{i}^{[l]})\right\|=\max_{j\in
1,\dots,(n_{S})}\frac{|R_{i}(S_{i}^{[l]})_{j}|}{\max(\text{abs},\,\text{rel}(S_{i}^{[l]})_{j})}<\tau,$
(8)

where $n_{S}$ is the dimension of $S_{i}$ and we choose $\tau=0.1$.
$M_{k}\approx\frac{\partial R_{i}(S_{i}^{[l]})}{\partial S_{i}^{[l]}}$ is the
iteration matrix defined as

$M_{k}=\begin{bmatrix}I-h_{k}\gamma\frac{\partial f(x_{k},y_{k},u)}{\partial
x_{k}}&-h_{k}\gamma\frac{\partial f(x_{k},y_{k},u)}{\partial y_{k}}\\\
-\frac{\partial g(x_{k},y_{k},u)}{\partial x_{k}}&-\frac{\partial
g(x_{k},y_{k},u)}{\partial y_{k}}\end{bmatrix},$ (9)

where $I$ is the identity matrix. For fixed step size ESDIRK methods, we
update $M_{k}$ at every integration step and we use its LU factorizations to
solve for all the stages in (7). Since the ESDIRK integration schemes are
stiffly accurate Runge-Kutta methods, we directly obtain the next differential
and algebraic states as

$x_{k+1}=X_{n_{s}}^{v_{n_{s},k}-1},\quad y_{k+1}=Y_{n_{s}}^{v_{n_{s},k}-1}.$
(10)

We apply stage value predictors (SVPs) of the form

$S_{i}^{[0]}=\alpha_{i}(r_{k})s_{k-1}+\sum_{j=2}^{n_{s}}\beta(r_{k})_{ij}\hat{S_{j}},$
(11)

to provide initial guesses for the Newton-type iterations in (7).
$\alpha(r_{k})\in\mathbb{R}^{n_{s}-1}$ and
$\beta(r_{k})\in\mathbb{R}^{n_{s}-1\times n_{s}-1}$ are the predictor
coefficients computed using $r_{k}=h_{k}/h_{k-1}$.
$\hat{S_{j}}=\hat{S_{j}}^{[v_{j,k-1}-1]}=\begin{bmatrix}\hat{X_{j}}^{[v_{j,k-1}-1]};&\hat{Y_{j}}^{[v_{j,k-1}-1]}\end{bmatrix}$
for $j=2,\dots n_{s}$ represent the previously converged differential and
algebraic stages and $s_{k-1}=\begin{bmatrix}x_{k-1};&y_{k-1}\end{bmatrix}$.
We apply the order conditions in [12] to compute the predictor coefficients
$\alpha(r_{k})$ and $\beta(r_{k})$, and we set $r_{k}=1$, as we only consider
fixed step-size ESDIRK methods. [7] shows the predictor coefficients for
ESDIRK12 and ESDIRK23. For the first step, $k=0$, we apply the trivial
predictor $S^{[0]}_{i}=s_{k}=\begin{bmatrix}x_{k};&y_{k}\end{bmatrix}$.

###  III-B Iterated IND sensitivity generation

We differentiate the adaptively generated discretization scheme of the ESDIRK
methods. For fixed step-size ESDIRK methods, the adaptive components are the
LU factorizations of the iteration matrices, $M_{k}$, the number of Newton-
type iterations for all stages, $v_{i,k}$, and the sequence of Newton-type
iterates, $S_{i}^{[l]}$ for $l=0,\dots,v_{i,k}-1$ and for $i=2,\dots,n_{s}$.
We define the sensitivities of the combined differential and algebraic states
with respect to initial states and inputs as

$\frac{\partial s_{k}}{\partial s_{0}}=\begin{bmatrix}\frac{\partial
x_{k}}{\partial x_{0}}&\frac{\partial x_{k}}{\partial y_{0}}\\\ \frac{\partial
y_{k}}{\partial x_{0}}&\frac{\partial y_{k}}{\partial
y_{0}}\end{bmatrix},\quad\frac{\partial s_{k}}{\partial
u}=\begin{bmatrix}\frac{\partial x_{k}}{\partial u}\\\ \frac{\partial
y_{k}}{\partial u}\end{bmatrix},$ (12)

with initial conditions $\frac{\partial s_{0}}{\partial s_{0}}=I$ and
$\frac{\partial s_{0}}{\partial u}=0$. We express the implementation of the
ESDIRK integration schemes in terms of the elementary operations

$\displaystyle S_{i}^{[0]}$
$\displaystyle=\phi^{\mathrm{SVP}}_{i}(s_{k-1},\hat{S}),$ (13a) $\displaystyle
S_{i}^{[l+1]}$ $\displaystyle=\phi^{\mathrm{it}}_{i}(S_{i}^{[l]},\psi_{i},u)$
$\displaystyle=S_{i}^{[l]}-M_{k}^{-1}R_{i}(S_{i}^{[l]},\psi_{i},u),$ (13b)
$\displaystyle s_{k+1}$
$\displaystyle=\phi^{n_{s}}(S_{n_{s}}^{[v_{n_{s},k}-1]})=S_{n_{s}}^{[v_{n_{s},k}-1]},$
(13c)

where (13a) represents the SVPs for stage $i$ with
$\hat{S}~{}=~{}\begin{bmatrix}\hat{S}_{2}^{[v_{2,k-1}-1]};&\dots;&s_{k}\end{bmatrix}$.
Eq. (13b) is the Newton-type scheme in (7) and (13c) represent the (10).

####  III-B1 State sensitivity

We construct the state sensitivities by computing the partial derivatives of
(13) with respect to the initial differential and algebraic states, $s_{0}$.
The sensitivities related to (13a) are

$\frac{\partial S_{i}^{[0]}}{\partial s_{0}}=\alpha_{i}(r_{k})\frac{\partial
s_{k-1}}{\partial
s_{0}}+\sum_{j=2}^{n_{s}}\beta_{ij}(r_{k})\frac{\partial\hat{S}_{j}}{\partial
s_{0}},$ (14)

for $k\geq 1$ and $\frac{\partial S_{i}^{[0]}}{\partial s_{0}}=\frac{\partial
s_{k}}{\partial s_{0}}$ for $k=0$. For (13b), we determine the state
sensitivities using the Newton-type scheme

$\frac{\partial S_{i}^{[l+1]}}{\partial s_{0}}=\frac{\partial
S_{i}^{[l]}}{\partial s_{0}}-M_{k}^{-1}\frac{\partial
R_{i}(S_{i}^{[l]},\psi_{i},u)}{\partial s_{0}},$ (15)

where $\frac{\partial R_{i}(S_{i}^{[l]},\psi_{i},u)}{\partial s_{0}}$ is

$\frac{\partial R_{i}(S_{i}^{[l]},\psi_{i},u)}{\partial s_{0}}=\frac{\partial
R_{i}}{\partial S_{i}^{[l]}}\frac{\partial S_{i}^{[l]}}{\partial
s_{0}}+\frac{\partial R_{i}}{\partial\psi_{i}}\frac{\partial\psi_{i}}{\partial
s_{0}}.$ (16)

The terms in (16) are

$\displaystyle\frac{\partial R_{i}}{\partial S_{i}^{[l]}}$
$\displaystyle=\begin{bmatrix}I-h_{k}\gamma\frac{\partial
f(X_{i}^{[l]},Y_{i}^{[l]},u)}{\partial X_{i}^{[l]}}&-h_{k}\gamma\frac{\partial
f(X_{i}^{[l]},Y_{i}^{[l]},u)}{\partial Y_{i}^{[l]}}\\\ -\frac{\partial
g(X_{i}^{[l]},Y_{i}^{[l]},u)}{\partial X_{i}^{[l]}}&-\frac{\partial
g(X_{i}^{[l]},Y_{i}^{[l]},u)}{\partial Y_{i}^{[l]}}\end{bmatrix},$ (17a)
$\displaystyle\frac{\partial R_{i}}{\partial\psi_{i}}$
$\displaystyle=\begin{bmatrix}-I\\\
0\end{bmatrix},\quad\frac{\partial\psi_{i}}{\partial
s_{0}}=\begin{bmatrix}\frac{\partial\psi_{i}}{\partial
x_{0}}&\frac{\partial\psi_{i}}{\partial y_{0}}\end{bmatrix}.$ (17b)

The partial derivatives in $\frac{\partial\psi_{i}}{\partial s_{0}}$ uses the
Jacobians of the previously converged stages for iteration $k$ as

$\displaystyle\frac{\partial\psi_{i}}{\partial x_{0}}$
$\displaystyle=\frac{\partial x_{k}}{\partial
x_{0}}+h_{k}\sum_{j=1}^{i-1}a_{ij}\bigg{(}\frac{\partial
f(X_{j},Y_{j},u)}{\partial S_{j}}\frac{\partial S_{j}}{\partial
x_{0}}\bigg{)}$ (18a) $\displaystyle\frac{\partial\psi_{i}}{\partial y_{0}}$
$\displaystyle=\frac{\partial x_{k}}{\partial
y_{0}}+h_{k}\sum_{j=1}^{i-1}a_{ij}\bigg{(}\frac{\partial
f(X_{j},Y_{j},u)}{\partial S_{j}}\frac{\partial S_{j}}{\partial
y_{0}}\bigg{)},$ (18b)

where

$\frac{\partial f(X_{j},Y_{j},u)}{\partial
S_{j}}=\begin{bmatrix}\frac{\partial f(X_{j},Y_{j},u)}{\partial
X_{j}}&\frac{\partial f(X_{j},Y_{j},u)}{\partial Y_{j}}\end{bmatrix}.$ (19)

Finally, due to the ESDIRK methods being stiffly accurate, we directly obtain
the sensitivities for (13c) as

$\frac{\partial x_{k+1}}{\partial x_{0}}=\frac{\partial
X_{n_{s}}^{v_{n_{s}},k-1}}{\partial x_{0}},\quad\frac{\partial
x_{k+1}}{\partial y_{0}}=\frac{\partial X_{n_{s}}^{v_{n_{s}},k-1}}{\partial
y_{0}}.$ (20)

####  III-B2 Input sensitivity

Similarly to the state sensitivities, we generate the input sensitivities by
computing the partial derivatives of (13) with respect to the vector of
inputs. For the initial guess of the Newton-type scheme we derive the
sensitivities at $k\geq 1$ as

$\frac{\partial S_{i}^{[0]}}{\partial u}=\alpha_{i}(r_{k})\frac{\partial
s_{k-1}}{\partial
u}+\sum_{j=2}^{n_{s}}\beta_{ij}(r_{k})\frac{\partial\hat{S}_{j}}{\partial u},$
(21)

with $\frac{\partial S_{i}^{[0]}}{\partial u}=\frac{\partial s_{k}}{\partial
u}$ for $k=0$. For the iterations in (13b) we compute the input sensitivity

$\frac{\partial S_{i}^{[l+1]}}{\partial u}=\frac{\partial
S_{i}^{[l]}}{\partial u}-M_{k}^{-1}\frac{\partial
R_{i}(S_{i}^{[l]},\psi_{i},u)}{\partial u},$ (22)

where the partial derivative of the residual function is

$\frac{\partial R_{i}(S_{i}^{[l]},\psi_{i},u)}{\partial u}=\frac{\partial
R_{i}}{\partial S_{i}^{[l]}}\frac{\partial S_{i}^{[l]}}{\partial
u}+\frac{\partial R_{i}}{\partial\psi_{i}}\frac{\partial\psi_{i}}{\partial
u}+\frac{\partial R_{i}}{\partial u}.$ (23)

The expressions $\frac{\partial R_{i}}{\partial S_{i}^{[l]}}$ and
$\frac{\partial R_{i}}{\partial\psi_{i}}$ are shown in (17) while we calculate
the other terms as

$\displaystyle\frac{\partial\psi_{i}}{\partial u}$
$\displaystyle=\frac{\partial x_{k}}{\partial u}$
$\displaystyle+h_{k}\sum_{j=1}^{i-1}a_{ij}\bigg{(}\frac{\partial
f(X_{j},Y_{j},u)}{\partial S_{j}}\frac{\partial S_{j}}{\partial
u}+\frac{\partial f(X_{j},Y_{j},u)}{\partial u}\bigg{)},$ (24a)
$\displaystyle\frac{\partial R_{i}}{\partial u}$
$\displaystyle=-\begin{bmatrix}h_{k}\gamma\frac{\partial
f(X_{i}^{[l]},Y_{i}^{[l]},u)}{\partial u};&\frac{\partial
g(X_{i}^{[l]},Y_{i}^{[l]},u)}{\partial u}\end{bmatrix}.$ (24b)

We directly obtain the input sensitivities for the next state as
$\frac{\partial x_{k+1}}{\partial u}=\frac{\partial
X_{n_{s}}^{v_{n_{s}},k-1}}{\partial u}$.

###  III-C Implicit-explicit SDAE solver

We solve (1) using an implicit-explicit integration scheme [1]. The implicit-
explicit integration scheme applies $M$-steps in each interval
$t\in[t_{k},t_{k+1}]$ for $k=0,\dots,N-1$ and is described by

$\displaystyle t_{k,n+1}$ $\displaystyle=t_{k,n}+\Delta t_{k,n},$ (25a)
$\displaystyle x_{k,n+1}$
$\displaystyle=x_{k,n}+f(t_{k,n+1},x_{k,n+1},y_{k,n+1},u_{k},d_{k})\Delta
t_{k,n}$ $\displaystyle+\sigma\Delta\omega_{k,n},$ (25b) $\displaystyle 0$
$\displaystyle=g(t_{k,n+1},x_{k,n+1},y_{k,n+1},u_{k},d_{k}),$ (25c)

for $n=0,\dots,M-1$ and where $t_{k,0}=t_{k}$, $x_{k,0}=x_{k}$,
$y_{k,0}=y_{k}$, and $\Delta\omega_{k,n}\sim N_{iid}(0,I\Delta t_{k,n})$. The
integration step-size is $\Delta t_{k,n}=t_{k,n+1}-t_{k,n}$ and we obtain the
step at $t_{k+1}$ as

$t_{k+1}=t_{k,M},\quad x_{k+1}=x_{k,M},\quad y_{k+1}=y_{k,M}.$ (26)

In each iteration $n$ in (25), we solve the system of nonlinear equations,
$R(s_{k,n+1})=0$. We define $R(s_{k,n+1})$ as

$R(s_{k,n+1})=\begin{bmatrix}x_{k,n+1}-x_{k,n}-f_{k,n}\Delta
t_{k,n}-\sigma\omega_{k,n}\\\ g_{k,n+1}\end{bmatrix},$ (27)

where $s_{k,n+1}=\begin{bmatrix}x_{k,n+1};&y_{k,n+1}\end{bmatrix}$ and

$\displaystyle f_{k,n}$
$\displaystyle=f(t_{k,n+1},x_{k,n+1},y_{k,n+1},u_{k},d_{k}),$ (28a)
$\displaystyle g_{k,n+1}$
$\displaystyle=g(t_{k,n+1},x_{k,n+1},y_{k,n+1},u_{k},d_{k}).$ (28b)

We solve (27) using an exact Newton scheme

$\displaystyle\frac{\partial R(s_{k,n+1}^{[l]})}{\partial
s_{k,n+1}^{[l]}}\Delta s_{k,n+1}^{[l]}$ $\displaystyle=-R(s_{k,n+1}^{[l]}),$
(29a) $\displaystyle s_{k,n+1}^{[l+1]}$ $\displaystyle=s_{k,n+1}^{[l]}+\Delta
s_{k,n+1}^{[l]},$ (29b)

and use a similar convergence criteria as in (8).

##  IV Nonlinear model predictive control

This section presents the implementations of the CD-EKF and the OCP. We
combine the CD-EKF and OCP into an NMPC such that we compute one control
interval between each measurement. We assume the vector of disturbances is
piece-wise constant between these measurements. We, therefore, apply the zero-
order-hold parametrizations

$u(t)=u_{k},\quad d(t)=d_{k}\quad\text{for}\quad t\in[t_{k},\,t_{k+1}],$ (30)

where the sampling time is $T_{s}=t_{k+1}-t_{k}$. $u_{k}$ represents the
implemented control signal computed by the OCP at time $t_{k}$.

###  IV-A Continuous-discrete extended Kalman filter

We implement the CD-EKF using the ESDIRK integration scheme for DAEs as in
[9]. The CD-EKF computes the filtered differential states, $\hat{x}_{k|k}$,
algebraic states, $\hat{y}_{k|k}$, and differential state covariance,
$P_{k|k}$, at time $t_{k}$ based on the measurements, $y^{m}_{k}$, and the
previous implemented control, $u_{k-1}$. We use the filtered states as initial
conditions for the OCP. The implemented control signal $u_{k}$ is used to
compute the one-step prediction at time $t_{k+1}$ in the differential states,
$\hat{x}_{k+1|k}$, algebraic states, $\hat{y}_{k+1|k}$, and the differential
state covariance, $P_{k+1|k}$. The initial states and covariance are

$\hat{x}_{0|-1}=x_{0},\quad\hat{y}_{0|-1}=y_{0},\quad P_{0|-1}=P_{0}.$ (31)

####  IV-A1 Filtering

We compute the innovation, $e_{k}$, as

$e_{k}=y_{k}^{m}-\hat{y}^{m}_{k|k-1},$ (32)

using the estimated measurement at time $t_{k}$

$\hat{y}^{m}_{k|k-1}=m(t_{k},\hat{x}_{k|k-1},\hat{y}_{k|k-1},u_{k-1},d_{k}).$
(33)

The covariance of the innovation, $R_{e,k}$, is computed as

$\displaystyle R_{e,k}$ $\displaystyle=C_{k}P_{k|k-1}C_{k}^{\prime}+R_{k},$
(34a) $\displaystyle C_{k}$ $\displaystyle=\frac{\partial m_{k|k-1}}{\partial
x}+\frac{\partial m_{k|k-1}}{\partial
y}\frac{\partial\hat{y}_{k|k-1}}{\partial x}.$ (34b)

The sensitivity $\frac{\partial\hat{y}_{k|k-1}}{\partial x}$ is computed by
solving

$\frac{\partial g_{k|k-1}}{\partial y}\frac{\partial y_{k|k-1}}{\partial
x}=-\frac{\partial g_{k|k-1}}{\partial x}.$ (35)

In these expressions, we used the notation

$\displaystyle m_{k|k-1}$
$\displaystyle=m(t_{k},\hat{x}_{k|k-1},\hat{y}_{k|k-1},u_{k-1},d_{k}),$ (36a)
$\displaystyle g_{k|k-1}$
$\displaystyle=g(t_{k},\hat{x}_{k|k-1},\hat{y}_{k|k-1},u_{k-1},d_{k}).$ (36b)

We combine (32)-(34) to compute the filtered state

$\hat{x}_{k|k}=\hat{x}_{k|k-1}+K_{k}e_{k},$ (37)

and the filtered covariance as

$P_{k|k}=(I-K_{k}C_{k})P_{k|k-1}(I-K_{k}C_{k})^{\prime}+K_{k}R_{k}K_{k}^{\prime},$
(38)

with $K_{k}$ being the Kalman gain defined as

$K_{k}=P_{k|k-1}C_{k}^{\prime}R_{e,k}^{-1}.$ (39)

Finally, to compute the filtered algebraic state, $\hat{y}_{k|k}$, we solve
the system of nonlinear equations

$0=g(t_{k},\hat{x}_{k|k},\hat{y}_{k|k},u_{k-1},d_{k}),$ (40)

using an exact Newton scheme.

####  IV-A2 Prediction

We compute one-step predictions, $\hat{x}_{k+1|k}=\hat{x}_{k}(t_{k+1})$ and
$\hat{y}_{k+1|k}=\hat{y}_{k}(t_{k+1})$, by solving

$\displaystyle\frac{d\hat{x}_{k}(t)}{dt}$
$\displaystyle=f(\hat{x}_{k}(t),\hat{y}_{k}(t),u_{k},d_{k}),$ (41a)
$\displaystyle 0$
$\displaystyle=g(\hat{x}_{k}(t),\hat{y}_{k}(t),u_{k},d_{k}),$ (41b)

for $t\in[t_{k},\,t_{k+1}]$, with initial conditions
$\hat{x}_{k}(t)=\hat{x}_{k|k}$ and $\hat{y}_{k}(t)=\hat{y}_{k|k}$ using an
ESDIRK method. [9] suggest computing the predicted differential state
covariance, $P_{k+1|k}=P_{k}(t_{k+1})$, as the solution to

$\begin{split}P_{k}(t)&=\Phi_{xx}(t,t_{k})P_{k}(t_{k})\Phi_{xx}^{\prime}(t,t_{k})\\\
&+\int_{t_{k}}^{t}\Phi_{xx}(t,s)\sigma\sigma^{\prime}\Phi_{xx}(t,s)ds,\end{split}$
(42)

at time $t=t_{k+1}$, with the initial condition $P_{k}(t_{k})=P_{k|k}$.
$\Phi_{xx}(t,s)=\frac{\partial\hat{x}_{k}(t)}{\partial\hat{x}_{k}(s)}$ is the
differential state sensitivity which we compute using the iterated IND
approach.

###  IV-B Optimal control problem

We consider an OCP that penalizes the tracking error, $z(t)-\bar{z}(t)$, and
the input rate-of-movement, $\Delta u_{k}=u_{k}-u_{k-1}$, in a least-squares
formulation. We solve the OCP at time $t_{k}$ with prediction horizon $T_{N}$
using the filtered differential state, $\hat{x}_{k|k}$, as initial conditions.
We divide the control horizon, $[t_{k},\,t_{k}+T_{N}]$, into $N=T_{N}/T_{s}$
equally spaced subintervals, $[t_{k+j},\,t_{k+j+1}]$, with
$j\in\mathcal{N}=1,\dots,N-1$. The OCP is

$\displaystyle\min_{x,y,u}\quad$ $\displaystyle\phi_{k},$ (43a) $\displaystyle
s.t.$ $\displaystyle x(t_{k})=\hat{x}_{k|k},$ (43b)
$\displaystyle\frac{dx(t)}{dt}=f(t,x,y,u,d),$ $\displaystyle
t\in[t_{k},\,t_{k}+T_{N}],$ (43c) $\displaystyle 0=g(t,x,y,u,d),$
$\displaystyle t\in[t_{k},\,t_{k}+T_{N}],$ (43d) $\displaystyle
z(t)=h(t,x,y,u,d),$ (43e) $\displaystyle u(t)=u_{k+j},\quad j\in\mathcal{N},$
$\displaystyle t\in[t_{k+j},\,t_{k+j+1}],$ (43f) $\displaystyle
d(t)=d_{k+j},\quad j\in\mathcal{N},$ $\displaystyle
t\in[t_{k+j},\,t_{k+j+1}],$ (43g) $\displaystyle u_{\min}\leq u_{k+j}\leq
u_{\max},$ $\displaystyle j\in\mathcal{N},$ (43h)

where $\phi_{k}=\phi_{k}([x(t);y(t);u(t)]_{t_{k}}^{t_{k}+T_{N}})$ is the
objective, and $f(t,x,y,u,d)=f(t,x(t),y(t),u(t),d(t))$ and
$g(t,x,y,u,d)=g(t,x(t),y(t),u(t),d(t))$ are the DAEs. We formulate the
objective $\phi_{k}=\phi_{z,k}+\phi_{\Delta u,k}+\phi_{N,k}$ with the
objectives function terms

$\displaystyle\phi_{z,k}$
$\displaystyle=\frac{1}{2}\int_{t_{k}}^{t_{k}+T_{N}}\left\|z(t)-\bar{z}(t)\right\|_{Q_{z}}^{2}dt,$
(44a) $\displaystyle\phi_{\Delta u,k}$
$\displaystyle=\frac{1}{2}\sum_{j=0}^{N-1}\left\|\Delta
u_{k+j}\right\|_{\bar{Q}_{\Delta u}}^{2},$ (44b) $\displaystyle\phi_{N,k}$
$\displaystyle=\frac{1}{2}\left\|z(t_{k}+T_{N})-\bar{z}(t_{k}+T_{N})\right\|_{\bar{Q}_{z}}^{2}.$
(44c)

$\phi_{z,k}$ penalizes the tracking error $z(t)-\bar{z}(t)$, with $\bar{z}(t)$
being the setpoint and $Q_{z}$ being the weight matrix. $\phi_{\Delta u,k}$
penalizes the input rate-of-movement with weight $\bar{Q}_{\Delta u}=Q_{\Delta
u}/T_{s}$. Finally, $\phi_{N,k}$ penalizes the tracking error at the end of
the horizon with the scaled weight matrix $\bar{Q}_{z}=Q_{z}/T_{s}$.

We apply the direct multiple shooting discretization scheme in [13] to
transcribe (43)-(44) into a NLP of the form

$\displaystyle\min_{w}\quad$ $\displaystyle\phi,$ (45a) $\displaystyle s.t.$
$\displaystyle w^{x}_{0}=\hat{x}_{k|k},$ (45b) $\displaystyle
R(w^{x}_{j+1},w^{x}_{j},w^{y}_{j},u_{j},d_{j})=0,\quad$ $\displaystyle
j\in\mathcal{N},$ (45c) $\displaystyle u_{\min}\leq u_{k+j}\leq u_{\max},\;$
$\displaystyle j\in\mathcal{N},$ (45d)

where $R(\cdot)$ is the ESDIRK integration scheme and the decision variables
are

$w=(w^{x}_{0},\,w^{y}_{0},\,u_{0},\,w^{x}_{1},\,w^{y}_{1},\,u_{1},\,\dots,u_{N-1},\,w^{x}_{N}).$
(46)

The $j\in\mathcal{N}$ are relative to $t_{k}$ and we initialize $w^{y}_{0}$
using the filtered algebraic state, $w^{y}_{0}=\hat{y}_{k|k}$. Since (44a)
depends continuously on the states, we discretize it by formulating it as an
IVP, which we solve simultaneously with the ESDIRK integration scheme. To
allow temporary violation of the algebraic constraints in the ESDIRK
integration scheme, we solve the relaxed DAEs [13] on each subinterval,

$\displaystyle\frac{dx_{j}(t)}{dt}$
$\displaystyle=f(t,x_{j}(t),y_{j}(t),u_{j},d_{j}),$ (47a) $\displaystyle 0$
$\displaystyle=g(t,x_{j}(t),y_{j}(t),u_{j},d_{j})$ $\displaystyle-
p_{j}(t)g(t_{j},w^{x}_{j},w^{y}_{j},u_{j},d_{j}),$ (47b)

with initial conditions $x_{j}(t_{j})=w^{x}_{j}$ and $y_{j}(t_{j})=w^{y}_{j}$.
$p_{j}(t)$ is a scalar function such that $p_{j}(t_{j})=1$ and $p_{j}(t)>0$
for $t\in[t_{k+j},\,t_{k+j+1}]$ and we choose
$p_{j}(t)=\exp({-\eta\frac{t-t_{j}}{t_{j+1}-t_{j}}})$ with $\eta=1$.

We solve (45) using an SQP solver with a BFGS for updating the Lagrangian
Hessian approximation. We denote the equality constraints (45b)-(45c) as
$b(w)$, and the gradient, $\frac{\partial b(w)}{\partial w}$, required by the
SQP solver, is computed using the sensitivity information $\frac{\partial
R(w)}{\partial w}$. As a result, the SQP solver requires both state and input
sensitivities. We implement the first control signal at time $t_{k}$,
$u^{*}_{0}$, from the solution $w^{*}$ to the system. For the first iteration
of the NMPC, $k=0$, we initialize $w$ using the initial consistent conditions,
$w_{0}=(x_{0},y_{0},u_{0},x_{0},y_{0},u_{0},\dots,x_{0})$, but for the
subsequent executions of the NLP, we initialize $w$ by shifting the previously
converged solution,
$w_{0}=(w^{x,*}_{1},w^{y,*}_{1},u_{1}^{*},\dots,w^{x,*}_{N},w^{y,*}_{N-1},u^{*}_{N-1},w^{x,*}_{N})$.

##  V Numerical case study

We apply the NMPC algorithm to an alkaline electrolyzer stack model. The model
we use is a simplified version of the one presented by [14], similar to [15].
The system of SDAEs are

$\displaystyle dT$
$\displaystyle=\frac{1}{C_{\mathrm{p,el}}}\Big{(}f_{\mathrm{in}}c_{\mathrm{P,lye}}(T_{\mathrm{in}}-T)$
$\displaystyle+n_{c}(U_{\mathrm{cell}}-U_{\mathrm{tn}})I-A_{s}h_{c}(T-T_{\mathrm{amb}})\Big{)}dt,$
(48a) $\displaystyle dT_{\mathrm{in}}$ $\displaystyle=\sigma d\omega,$ (48b)
$\displaystyle 0$
$\displaystyle=U_{\mathrm{cell}}-(U_{\mathrm{rev}}+U_{\mathrm{ohm}}(T,I)+U_{\mathrm{act}}(T,I)),$
(48c) $\displaystyle 0$
$\displaystyle=P_{\mathrm{in}}-n_{c}U_{\mathrm{cell}}I,$ (48d)

where

$\displaystyle U_{\mathrm{ohm}}(T,I)$
$\displaystyle=(r_{1}+r_{2}T)\frac{I}{A},$ (49a) $\displaystyle
U_{\mathrm{act}}(T,I)$
$\displaystyle=s\log\left(\left(t_{1}+\frac{t_{2}}{T}+\frac{t_{3}}{T^{2}}\right)\frac{I}{A}+1\right).$
(49b)

The SDAEs model the energy balance and electrochemical dynamics in the
electrolyzer stack. $x(t)=\begin{bmatrix}T;&T_{in}\end{bmatrix}$ being the
vector of differential variables representing the temperature in the
electrolyzer stack and the inlet temperature, respectively. In this system,
the inlet temperature is modeled as a stochastic differential equation with a
zero drift and where $\sigma=0.03$.
$y(t)=\begin{bmatrix}U_{\mathrm{cell}};&I\end{bmatrix}$ is the vector of
algebraic equations with $U_{\mathrm{cell}}$ being the cell voltage and $I$
being the current. $C_{\mathrm{p,el}}$, $c_{\mathrm{P,lye}}$,
$n_{c},U_{\mathrm{tn}}$, $A_{s}$, $h_{c}$, $U_{\mathrm{rev}}$, and $A$ are
parameters of the model and represent the overall stack heat capacity, lye
circulation heat capacity, number of electrolytic cells, thermoneutral
voltage, active transfer area, heat transfer coefficient, reversible voltage,
and electrodes area, respectively. $r_{1}$, $r_{2}$, $t_{1}$, $t_{2}$,
$t_{3}$, and $s$ are the parameters related to the electrochemical model. We
use the values for the parameters provided by [15]. The lye inlet flow rate
$u(t)=f_{in}$, is the input. The disturbances,
$d(t)=\begin{bmatrix}T_{\mathrm{amb}};&P_{\mathrm{in}}\end{bmatrix}$, are the
ambient temperature and the inlet power, respectively, with the values
$T_{\mathrm{amb}}=25$ $C^{\circ}$ and $P_{\mathrm{in}}=2$ MW. We assume these
disturbances to be known by the NMPC.

We measure the temperature, i.e., $y_{k}^{m}=T_{k}+v_{k}$ in (1c), where
$v_{k}$ is the measurement noise with covariance $R_{k}=1$. We sample the
measurements at every $T_{s}=4$ min. The NMPC tracks a time-varying setpoint
for the temperature, $T$. Hence, The output equation in (1d) is $z(t)=T$. For
the OCP, we apply the control horizon $T_{N}=100$ min using $N=25$ steps. The
weight matrices for tracking error and input-rate-of-movement penalties are
$Q_{z}=10$ and $Q_{\Delta u}=0.1$, respectively. We constraint the inlet flow
to $2\;\text{kg/s}\leq f_{in}\leq 10\;\text{kg/s}$ and apply the ESDIRK34 with
fixed step size $h_{k}=h=0.2T_{s}$ in both the OCP and CD-EKF. We use the
implicit-explicit method to perform stochastic simulations of the closed-loop
system.

Fig 1 shows the closed-loop simulation of the alkaline electrolyzer stack,
tracking a time-varying setpoint $\bar{z}(t)$ varying between 75∘C and 60∘C.
The result shows that the NMPC effectively follows the setpoint, even when the
initial estimate of the inlet temperature, $\hat{T}_{\mathrm{in}}$, deviates
from the actual state. Improved tracking accuracy is observed as the
estimation error in $T_{\mathrm{in}}$ diminishes.

Figure 1: Closed-loop simulation of the alkaline electrolyzer model with NMPC
based on ESDIRK34 integration.

##  VI Conclusions

We present an NMPC algorithm for setpoint tracking in systems described by
continuous-discrete semi-explicit SDAEs of index 1. The NMPC consists of a CD-
EKF and an OCP discretized by direct multiple shooting. We apply ESDIRK
integration schemes to solve the IVPs in the CD-EKF and OCP. The ESDIRK
implementation uses an iterative IND approach to compute sensitivities. These
sensitivities are used to provide accurate gradient information to the OCPs
and to compute the covariance in the CD-EKF. We apply the NMPC to an alkaline
electrolyzer stack model to track a time-varying setpoint for temperature.

## References

  * [1] T. K. S. Ritschel and J. B. Jørgensen, “Nonlinear Model Predictive Control for Disturbance Rejection in Isoenergetic-isochoric Flash Processes,” in _IFACPapersOnLine_ , vol. 52, no. 1, 2019, pp. 796–801.
  * [2] L. T. Biegler, S. L. Campbell, and V. Mehrmann, _Control and Optimization with Differential-Algebraic Constraints_. SIAM: Society for Industrial and Applied Mathematics, 2012.
  * [3] U. M. Ascher and L. R. Petzold, _Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations_. SIAM: Society for Industrial and Applied Mathematics, Jan. 1998.
  * [4] M. R. Kristensen, J. B. Jørgensen, P. G. Thomsen, and S. B. Jørgensen, “An ESDIRK method with sensitivity analysis capabilities,” _Computers & Chemical Engineering_, vol. 28, no. 12, pp. 2695–2707, 2004.
  * [5] H. G. Bock, “Numerical Treatment of Inverse Problems in Chemical Reaction Kinetics,” in _Modelling of Chemical Reaction Systems_ , K. H. Ebert, P. Deuflhard, and W. Jäger, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 1981, pp. 102–125.
  * [6] J. Albersmeyer, “Adjoint based algorithms and numerical methods for sensitivity generation and optimal control of large scale dynamic systems,” Ph.D. dissertation, Ruprecht-Karls-Universität Heidelberg, Dec. 2010.
  * [7] A. H. D. Christensen and J. B. Jørgensen, “A Comparative Study of Sensitivity Computations in ESDIRK-Based Optimal Control Problems,” in _Accepted for the 22nd European Control Conference (ECC24)_ , 2024.
  * [8] J. B. Jørgensen, “A critical discussion of the continuous-discrete extended Kalman filter,” in _European Congress of Chemical Engineering-6_ , Copenhagen, Denmark, 2007.
  * [9] J. B. Jørgensen, M. R. Kristensen, P. G. Thomsen, and H. Madsen, _Assessment and Future Directions of Nonlinear Model Predictive Control_ , ser. Lecture Notes in Control and Information Sciences. Springer Berlin, Heidelberg, 2007, vol. 358, ch. New extended Kalman filter algorithms for stochastic differential algebraic equations, pp. 359–366.
  * [10] M. R. Kristensen, J. B. Jørgensen, P. G. Thomsen, M. L. Michelsen, and S. B. Jørgensen, “Sensitivity analysis in index-1 differential algebraic equations by ESDIRK methods,” in _16th IFAC World Congress_ , vol. 38, no. 1, Prague, Czech Republic, 2005, pp. 212–217, 16th IFAC World Congress.
  * [11] J. B. Jørgensen, M. R. Kristensen, and P. G. Thomsen, “A Family of ESDIRK Integration Methods,” _arXiv preprint arXiv:1803.01613_ , 2018.
  * [12] I. Higueras and T. Roldán, “Starting Algorithms for a Class of RK Methods for Index-2 DAEs,” _Computers & Mathematics with Applications_, vol. 49, pp. 1081–1099, 2005.
  * [13] H. G. Bock, M. Diehl, P. Kühl, E. Kostina, J. P. Schlöder, and L. Wirsching, _Assessment and Future Directions of Nonlinear Model Predictive Control_ , ser. Lecture Notes in Control and Information Sciences. Springer Berlin, Heidelberg, 2007, vol. 358, ch. Numerical Methods for Efficient and Fast Nonlinear Model Predictive Control, pp. 163–179.
  * [14] N. Cantisani, J. Dovits, and J. B. Jørgensen, “Dynamic modeling of an alkaline electrolyzer plant for process simulation and optimization,” _arXiv preprint arXiv:2311.09882_ , 2023.
  * [15] M. Rizwan, V. Alstad, and J. Jäschke, “Design considerations for industrial water electrolyzer plants,” _International Journal of Hydrogen Energy_ , vol. 46, pp. 37 120–37 136, 2021.