# ALMol: Aligned Language-Molecule Translation LLMs through Offline Preference
Contrastive Optimisation

Dimitris Gkoumas
Queen Mary University of London, London, UK
<EMAIL_ADDRESS>

###### Abstract

The field of chemistry and Artificial Intelligence (AI) intersection is an
area of active research that aims to accelerate scientific discovery. The
integration of large language models (LLMs) with scientific modalities has
shown significant promise in this endeavour. However, challenges persist in
effectively addressing training efficacy and the out-of-distribution problem,
particularly as existing approaches rely on larger models and datasets. In
this context, we focus on machine language-molecule translation and deploy a
novel training approach called contrastive preference optimisation, which
avoids generating translations that are merely adequate but not perfect. To
ensure generalisability and mitigate memorisation effects, we conduct
experiments using only 10% of the data. Our results demonstrate that our
models achieve up to a 32% improvement compared to counterpart models. We also
introduce a scalable fine-grained evaluation methodology that accommodates
responsibility.

ALMol: Aligned Language-Molecule Translation LLMs through Offline Preference
Contrastive Optimisation

Dimitris Gkoumas Queen Mary University of London, London, UK
<EMAIL_ADDRESS>

##  1 Introduction

The world is facing unprecedented complexity in the form of global challenges
such as climate change, healthcare, and pandemics. Innovative scientific
solutions are urgently needed to address these challenges. Chemistry has been
at the forefront of developing such solutions, pioneering new drugs Ferguson
and Gray (2018), creating advanced materials Kippelen and Brédas (2009), or
enhancing chemical processes Zhong et al. (2023). However, these frontiers are
vast and require the involvement of Artificial Intelligence (AI) technology to
navigate them effectively.

Large language models (LLMs) have shown promising potential for accelerating
scientific discovery across various domains, including chemistry, biology, and
materials science Zhang et al. (2023); AI4Science and Quantum (2023).
Researchers have explored successful paradigms from natural language
processing (NLP) and multimodal representation learning to the chemistry
domain. One prevalent strategy involves extrapolating the inherent three-
dimensional structure of molecules to SMILES, representing them as character
strings, and subsequently learning language-molecule representations either in
separate yet coordinated spaces Edwards et al. (2022, 2021); Liu et al.
(2023a), in a joint space Liu et al. (2023b), or through hybrid approaches Luo
et al. (2023); Christofidellis et al. (2023). In light of the recent
significant advancements in the field, none of the above approaches
effectively tackle the inherent challenges in training such models. Instead,
they rely on sparse or noisy synthetic data, often necessitating exponentially
more data than is typically used in NLP tasks Edwards et al. (2024).

However, training on larger models and datasets does not necessarily guarantee
higher performance. A successful paradigm that augments the capabilities of
LLMs across multiple NLP tasks is Reinforcement Learning with Human Feedback
(RLHF) Ouyang et al. (2022). Although initially challenged by issues of
slowness and instability, recent research has addressed many of these
challenges by shifting towards closed-form losses that operate directly on
offline preferences Rafailov et al. (2024). RLHF has demonstrated superior
performance compared to standard minimising cross-entropy optimisation
approaches.

In this context, we address challenges related to effectively training robust
language models when integrated with scientific modalities. We deploy a novel
way of training LLMs for language-molecule translation that avoids generating
translations that are only adequate but not perfect, called contrastive
preference optimisation (CTO) Xu et al. (2024). CTO is based on offline
preferences instead of supervised fine-tuning, mimicking reference
translations. To ensure that our models can effectively generalise instead of
memorising patterns, we conduct experiments using only 10% of the L+M-24
dataset Edwards et al. (2024). Our contributions have as follows:

  * •

Our models achieve significant improvement in performance across a variety of
evaluation metrics compared with models trained on extensive in and out-of-
distribution data ($\S~{}\ref{sec:exp_results}$).

  * •

We demonstrate their robustness by experimenting with pivot versus minor
modalities ($\S~{}\ref{sec:exp_results}$).

  * •

We introduce a fine-grained evaluation methodology that is domain-independent
and enhances explainability ($\S~{}\ref{sec:eval_method}$).

##  2 Background

Reinforcement Learning with Human Feedback (RLHF) optimisation Ouyang et al.
(2022) operates with a triple dataset $\mathcal{D}=\\{x,y_{w},y_{l}\\}$, where
$y_{w}$ and $y_{l}$ represent preferred and dis-preferred outputs,
corresponding to input $x$, such that $y_{w}\succ y_{l}$ for $x$. The
probability of $y_{w}$ over $y_{l}$ in pairwise comparisons is typically
computed using the Bradley-Terry model Bradley and Terry (1952):

$p^{*}(y_{w}\succ y_{l}|x)=\sigma(r^{*}(x,y_{w})-r^{*}(x,y_{l}))$ (1)

where $\sigma$ is the logistic function, and $r^{*}$ denotes the reward
function that underlies the preferences.

As obtaining the reward directly from a human would be prohibitively
expensive, a reward model $r_{\phi}$ is trained to act as a surrogate by
minimising the negative log-likelihood of the preference data:

$\mathcal{L}(r_{\phi})=-\mathbb{E}_{(x,y_{w},y_{l})\sim\mathcal{D}}[\log\sigma(r_{\phi}(x,y_{w})-r_{\phi}(x,y_{l}))]$
(2)

Additionally, the Kullback-Leibler (KL) divergence between the outputs
generated by a reference model $\pi_{\text{ref}}$ and the parameterised
$\pi_{\theta}$ models serves as an additional reward signal, ensuring that the
generated responses closely align with the reference model. Consequently, an
optimal model $\pi_{\theta}$ is one that maximises:

$\mathbb{E}_{(x\in\mathcal{D},y\in\pi_{\theta})}[r_{\phi}(x,y)]-\beta\mathcal{D}_{\text{KL}}\bigr{(}\pi_{\theta}(y|x)||\pi_{\text{ref}}(y|x)\bigr{)}$
(3)

where $\beta$ is the temperature parameter typically $\in[0.1,0.5]$.

RLHF can present challenges due to its inherent slowness and instability,
especially in distributed settings Zheng et al. (2024). Recent work has
shifted towards closed-form losses to align LLMs with human preferences. Here,
we experiment with contrastive preference optimisation that adopts a closed-
form loss for RLHF.

##  3 Methodology

###  3.1 Task Formulation

Let $L$ denote the source language space and $M$ denote the target molecule
space. Consider a dataset $\mathcal{D}$, comprising pairs of source language
sequences $x$ represented in space $L$, and their corresponding target
molecule sequences $y$ in space $M$, represented as
$\mathcal{D}=\\{x^{(i)},y^{(i)}\\}_{i=1}^{N}$, where $N$ is the total number
of pairs. Here, we cast the problem on a machine translation task, aiming to
learn an optimal function $f:L\leftrightarrow M$ through a model
$\pi_{\theta}$ parameterised by $\theta$. We unify language and molecule
representations through instructional modelling to direct the translation
process, allowing for both $L\rightarrow M$ and $M\rightarrow L$ (please refer
to Appx. A).

###  3.2 Contrastive Preference Optimisation

Contrastive preference optimisation (CTO) Xu et al. (2024) addresses
challenges stemming from the inherent limitation in RLHF, as discussed in $\S$
2, and from the necessity of high-quality data. CTO is a general approximation
of Eq. 3 using a uniform reference model, which assumes equal likelihood for
all possible generated outputs:

$\displaystyle\mathcal{L}(\pi_{\theta};U)=-\mathbb{E}_{(x,y_{w},y_{l})\sim D}$
$\displaystyle\Biggl{[}\log\sigma\Bigl{(}\beta\log\pi_{\theta}(y_{w}|x)-\beta\log\pi_{\theta}(y_{l}|x)\Bigr{)}\Biggr{]}$
(4)

where $\pi_{\theta}$ is parameterised model by $\theta$ and $\beta$ hyper-
parameter (please refer $\S$ 2). Eq. 4 implies that the loss is calculated
based on how well the generated translations match this uniform distribution
of possible translations, rather than being biased towards any particular
translation. To maintain $\pi_{\theta}$ close to the preferred data
distribution, a behaviour cloning (BC) Hejna et al. (2023) regulariser is
introduced:

$\displaystyle\min_{\theta}\mathcal{L}(\pi_{\theta},U)\quad\text{s.t.}\quad$
$\displaystyle\mathbb{E}_{(x,y_{w})\sim
D}\Bigr{[}\mathbb{KL}\bigl{(}\pi_{w}(y_{w}|x)||\pi_{\theta}(y_{w}|x)\bigl{)}\Bigr{]}<\epsilon,$
(5)

Here, $\epsilon$ denotes a small positive constant, and $\mathbb{KL}$
signifies the Kullback-Leibler divergence. The regulariser is enhanced with an
additional SFT term on the preferred data, bolstering the CPO loss as:

$\mathcal{L}_{\text{CPO}}=\min_{\theta}\underbrace{\mathcal{L}(\pi_{\theta},U)}_{\mathcal{L}_{\text{prefer}}}-\underbrace{\mathbb{E}_{(x,y_{w})\sim
D}\bigr{[}\log\pi_{\theta}(y_{w}|x)\bigr{]}}_{\mathcal{L}_{\text{NLL}}}$ (6)

##  4 Experiments

###  4.1 Data

We conduct experiments on the L+M-24 benchmark dataset, which encompasses both
molecule and linguistic modalities Edwards et al. (2024). It is divided into
four categories, each with significant applications in small-molecule domain;
biomedical; light and electricity; human interaction and organoleptics; and
agriculture and industry. The training and validation subsets consist of
approximately 127k and 34k language-molecule pairs, respectively. Here, we
utilise 10% of these subsets for training and validation. To operationalise
CTO, we recreate a triples dataset consisting of preferred and dis-preferred
outputs (see $\S$ 2), where the former are the golden references and the
latter are generated from MolT5 Edwards et al. (2022). For testing, we sample
3k unseen pairs from the test subset of additional training data released for
the shared task at Language + Molecules $@$ ACL2024.

###  4.2 Bechmark Models

We compare our results with established language-molecule models as captured
in the literature:

  * •

TxtChem-T5 Christofidellis et al. (2023): A T5 model trained on both
linguistic and molecule modalities with a multi-task objective across various
datasets, including the CheBI-20 dataset Edwards et al. (2022), akin to
L+M-24.

  * •

Chem-LLM Zhang et al. (2024): An InternLM2-Base-7B model, trained on an
extensive chemical domain knowledge dataset, with the direct preference
optimisation objective Rafailov et al. (2024), achieves results comparable to
GPT-4.

  * •

Meditron Chen et al. (2023): A Meditron-7B model fine-tuned on the entire
L+M-24 for unidirectional language-molecule translation.

  * •

SFT-Meditron: We fine-tune Meditron-7B on a 10% subset of L+M-24 for bi-
directional machine language-molecule translation.

###  4.3 Experimental Settings

Here, we train Meditron with CTO on a 10% subset of L+M-24. We experiment with
both language and molecule weight initialisation obtained from Meditron
trained on the entire data Edwards et al. (2024). We refer to them as CTO-
Meditron${}_{\overrightarrow{Lan.}}$ and CTO-
Meditron${}_{\overrightarrow{Mol.}}$, respectively. We train the models with
QLoRA Dettmers et al. (2024) (more details in Appx C). For evaluation, we
adopt established metrics in Edwards et al. (2022).

###  4.4 Experiment Results

Table 2 presents a summary of the molecule-to-language results. We observed a
significant decrease in performance for benchmark models trained on extensive
data with SFT when tested on out-of-distribution data. Among the baseline
models, Meditron demonstrated the highest performance, likely due to its
training on the entire L+M-24 dataset utilised in our experiments. Training
Meditron with SFT for bi-directional language-molecule translation has
demonstrated neither effectiveness (see Table 1) nor efficiency (refer to
Appx. B). This suggests that the performance in our experiments is not
dependent on memorised patterns from Meditron trained on the entire dataset.
In contrast, our models trained with the CTO objective on only 10% of L+M-24
achieved a remarkable improvement in performance across diverse evaluation
metrics, up to 32% compared to Meditron trained on the entire dataset. This
improvement is consistent, as our model consistently enhances performance when
initialised from modality-agnostic weights, i.e., CTO-
Meditron${}_{\overrightarrow{Lan}}$ in Table 1.

Model | $\textbf{Blue-2}\uparrow$ | $\textbf{Blue-4}\uparrow$ | $\textbf{Rouge-1}\uparrow$ | $\textbf{Rouge-2}\uparrow$ | $\textbf{Rouge-L}\uparrow$ | $\textbf{METEOR}\uparrow$
---|---|---|---|---|---|---
TxtChem-T5 | 0.08 | 0.09 | 0.19 | 0.06 | 0.17 | 0.16
Chem-LLM | 0.03 | 0.00 | 0.11 | 0.02 | 0.09 | 0.14
Meditron | 0.42 | 0.30 | 0.63 | 0.47 | 0.49 | 0.54
SFT-Meditron | 0.37 | 0.26 | 0.54 | 0.39 | 0.38 | 0.60
CTO-Meditron${}_{\overrightarrow{Lan}}$ | 0.62 (+0.20) | 0.45 (+0.15) | 0.67 (+0.03) | 0.50 (+0.03) | 0.48 (-0.01) | 0.62 (+0.08)
CTO-Meditron${}_{\overrightarrow{Mol}}$ | 0.74 (+0.32) | 0.53 (+0.23) | 0.76 (+0.10) | 0.56 (+0.09) | 0.53 (+0.04) | 0.71(+0.17)

Table 1: Molecule-to-language translation results. Arrows next to metrics
indicate the higher value the better performance. Numbers in parentheses show
deviations from Meditron trained on the entire dataset.

We observed similar performance patterns for language-to-molecule translation
as reported in Table 2. However, even though our model achieved better
performance compared to Meditron when initialised with modality-agnostic
weights, it struggled to learn molecular patterns (see CTO-
Meditron${}_{\overrightarrow{Mol.}}$ in Table 2). This suggests that language
plays a pivotal role in the molecule modality. In the future, we aim to
explore more advanced initialised methods to address this challenge.

Model | $\textbf{BLEU}\uparrow$ | $\textbf{Exact}\uparrow$ | $\textbf{Levenshtein}\downarrow$ | $\textbf{MACCS FTS}\uparrow$ | $\textbf{RDK FTS}\uparrow$ | $\textbf{Morgan FTS}\uparrow$ | $\textbf{FCD}\downarrow$ | $\textbf{Validity}\uparrow$
---|---|---|---|---|---|---|---|---
TxtChem-T5 | 0.18 | 0.00 | 133.29 | 0.21 | 0.10 | 0.03 | 37.67 | 0.58
Chem-LLM | 0.04 | 0.00 | 732.74 | 0.00 | 0.00 | 0.00 | 59.44 | 0.19
Meditron | 0.43 | 0.00 | 66.16 | 0.35 | 0.29 | 0.19 | 13.64 | 0.57
SFT-Meditron | 0.30 | 0.00 | 186.99 | 0.70 | 0.62 | 11.14 0.41 |  | 0.98
CTO-Meditron${}_{\overrightarrow{Lan.}}$ | 0.71 (+0.28) | 0.00 | 42.65 (-23.51) | 0.78 (+0.43) | 0.70 (+0.41) | 0.48 (+0.29) | 4.19 (-9.45) | 1.00 (+0.43)
CTO-Meditron${}_{\overrightarrow{Mol.}}$ | 0.52 (+0.09) | 0.00 | 76.95 (+10.43) | 0.52 (+0.17) | 0.49 (+0.20) | 0.37 (+0.18) | 27.39 (+13.75) | 0.58 (+0.01)

Table 2: Language-to-molecule translation results. Arrows next to metrics
indicate whether higher or lower values denote better performance. Numbers in
parentheses show deviations from Meditron trained on the entire dataset.

###  4.5 Proposed Evaluation Methodology

Prior studies have utilised embedding representations, for assessing the
semantics in chemical-domain models Jaeger et al. (2018); Edwards et al.
(2021); Christofidellis et al. (2023). However, these approaches require
domain adaptation for out-of-distribution data Edwards et al. (2024) and might
lead to opaque and arbitrary outcomes Steck et al. (2024). We address these
limitations by introducing a scalable fine-grained evaluation methodology.

For molecule-to-language translation, we utilise scores from a natural
language inference (NLI) model111https://huggingface.co/cross-encoder/nli-
deberta-v3-base based on He et al. (2021).. Specifically, we obtain the
probabilities of predictions being entailed in references. Additionally, we
compute the length deviation of predictions from the references. Evaluation
results, as illustrated in Fig. 1, indicated that CTO-
Meditron${}_{\overrightarrow{Mol.}}$ excels at generating outputs with lengths
similar to the actual ones, whereas CTO-Meditron${}_{\overrightarrow{Lan.}}$
is most adept at capturing language semantics.

Figure 1: LLM hallucination for molecule-to-language translation. The upper
figure displays distributions of length deviations and the lower one
entailment probabilities in prediction-reference pairs. Length closer to zero
and higher probability indicate superior performance. Figure 2: LLM
hallucination for language-to-molecule translation. The upper figure displays
distributions of length deviations and the lower one CharacTER scores in
prediction-reference pairs. Lower CharacTER scores and length deviations
closer to zero indicate better performance.

For language-to-molecule translation, we employ a character-level metric,
CharacTER Wang et al. (2016), which measures the minimum number of character
edits needed to adjust a prediction until it matches the reference. Similarly,
we calculate length deviations in prediction-inference pairs. The evaluation
results in Fig. 2 show that CTO-Meditron${}_{\overrightarrow{Lan.}}$ is the
most reliable model.

##  5 Conclusion

This work address training efficacy and the out-of-distribution problem for
language models integrated with scientific modalities. We train models for
machine language-molecule translation on only 10% of available data with
contrastive preference optimisation which avoids generating translations that
are merely adequate but not perfect. We achieve significant improvement in
performance when compared with models trained on extensive in and out-of-the-
distribution data. Finally, we introduce a fine-grained evaluation methodology
that is domain-independent and aids explainability.

## References

  * AI4Science and Quantum (2023) Microsoft Research AI4Science and Microsoft Azure Quantum. 2023.  The impact of large language models on scientific discovery: a preliminary study using gpt-4.  _arXiv preprint arXiv:2311.07361_.
  * Bradley and Terry (1952) Ralph Allan Bradley and Milton E Terry. 1952.  Rank analysis of incomplete block designs: I. the method of paired comparisons.  _Biometrika_ , 39(3/4):324–345.
  * Chen et al. (2023) Zeming Chen, Alejandro Hernández Cano, Angelika Romanou, Antoine Bonnet, Kyle Matoba, Francesco Salvi, Matteo Pagliardini, Simin Fan, Andreas Köpf, Amirkeivan Mohtashami, et al. 2023.  Meditron-70b: Scaling medical pretraining for large language models.  _arXiv preprint arXiv:2311.16079_.
  * Christofidellis et al. (2023) Dimitrios Christofidellis, Giorgio Giannone, Jannis Born, Ole Winther, Teodoro Laino, and Matteo Manica. 2023.  Unifying molecular and textual representations via multi-task language modelling.  In _International Conference on Machine Learning_ , pages 6140–6157. PMLR.
  * Dettmers et al. (2024) Tim Dettmers, Artidoro Pagnoni, Ari Holtzman, and Luke Zettlemoyer. 2024.  Qlora: Efficient finetuning of quantized llms.  _Advances in Neural Information Processing Systems_ , 36.
  * Edwards et al. (2022) Carl Edwards, Tuan Lai, Kevin Ros, Garrett Honke, Kyunghyun Cho, and Heng Ji. 2022.  Translation between molecules and natural language.  In _Proceedings of the 2022 Conference on Empirical Methods in Natural Language Processing_ , pages 375–413.
  * Edwards et al. (2024) Carl Edwards, Qingyun Wang, Lawrence Zhao, and Heng Ji. 2024.  L+ m-24: Building a dataset for language+ molecules@ acl 2024.
  * Edwards et al. (2021) Carl Edwards, ChengXiang Zhai, and Heng Ji. 2021.  Text2mol: Cross-modal molecule retrieval with natural language queries.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 595–607.
  * Ferguson and Gray (2018) Fleur M Ferguson and Nathanael S Gray. 2018.  Kinase inhibitors: the road ahead.  _Nature reviews Drug discovery_ , 17(5):353–377.
  * He et al. (2021) Pengcheng He, Jianfeng Gao, and Weizhu Chen. 2021.  Debertav3: Improving deberta using electra-style pre-training with gradient-disentangled embedding sharing.  _arXiv preprint arXiv:2111.09543_.
  * Hejna et al. (2023) Joey Hejna, Rafael Rafailov, Harshit Sikchi, Chelsea Finn, Scott Niekum, W Bradley Knox, and Dorsa Sadigh. 2023.  Contrastive prefence learning: Learning from human feedback without rl.  _arXiv preprint arXiv:2310.13639_.
  * Jaeger et al. (2018) Sabrina Jaeger, Simone Fulle, and Samo Turk. 2018.  Mol2vec: unsupervised machine learning approach with chemical intuition.  _Journal of chemical information and modeling_ , 58(1):27–35.
  * Kippelen and Brédas (2009) Bernard Kippelen and Jean-Luc Brédas. 2009.  Organic photovoltaics.  _Energy & Environmental Science_, 2(3):251–261.
  * Liu et al. (2023a) Shengchao Liu, Weili Nie, Chengpeng Wang, Jiarui Lu, Zhuoran Qiao, Ling Liu, Jian Tang, Chaowei Xiao, and Animashree Anandkumar. 2023a.  Multi-modal molecule structure–text model for text-based retrieval and editing.  _Nature Machine Intelligence_ , 5(12):1447–1457.
  * Liu et al. (2023b) Zequn Liu, Wei Zhang, Yingce Xia, Lijun Wu, Shufang Xie, Tao Qin, Ming Zhang, and Tie-Yan Liu. 2023b.  Molxpt: Wrapping molecules with text for generative pre-training.  _arXiv preprint arXiv:2305.10688_.
  * Luo et al. (2023) Yizhen Luo, Kai Yang, Massimo Hong, Xingyi Liu, and Zaiqing Nie. 2023.  Molfm: A multimodal molecular foundation model.  _arXiv preprint arXiv:2307.09484_.
  * Ouyang et al. (2022) Long Ouyang, Jeffrey Wu, Xu Jiang, Diogo Almeida, Carroll Wainwright, Pamela Mishkin, Chong Zhang, Sandhini Agarwal, Katarina Slama, Alex Ray, et al. 2022.  Training language models to follow instructions with human feedback.  _Advances in neural information processing systems_ , 35:27730–27744.
  * Rafailov et al. (2024) Rafael Rafailov, Archit Sharma, Eric Mitchell, Christopher D Manning, Stefano Ermon, and Chelsea Finn. 2024.  Direct preference optimization: Your language model is secretly a reward model.  _Advances in Neural Information Processing Systems_ , 36.
  * Steck et al. (2024) Harald Steck, Chaitanya Ekanadham, and Nathan Kallus. 2024.  Is cosine-similarity of embeddings really about similarity?  _arXiv preprint arXiv:2403.05440_.
  * Wang et al. (2016) Weiyue Wang, Jan-Thorsten Peter, Hendrik Rosendahl, and Hermann Ney. 2016.  Character: Translation edit rate on character level.  In _Proceedings of the First Conference on Machine Translation: Volume 2, Shared Task Papers_ , pages 505–510.
  * Xu et al. (2024) Haoran Xu, Amr Sharaf, Yunmo Chen, Weiting Tan, Lingfeng Shen, Benjamin Van Durme, Kenton Murray, and Young Jin Kim. 2024.  Contrastive preference optimization: Pushing the boundaries of llm performance in machine translation.  _arXiv preprint arXiv:2401.08417_.
  * Zhang et al. (2024) Di Zhang, Wei Liu, Qian Tan, Jingdan Chen, Hang Yan, Yuliang Yan, Jiatong Li, Weiran Huang, Xiangyu Yue, Dongzhan Zhou, et al. 2024.  Chemllm: A chemical large language model.  _arXiv preprint arXiv:2402.06852_.
  * Zhang et al. (2023) X Zhang, L Wang, J Helwig, Y Luo, C Fu, Y Xie, M Liu, Y Lin, Z Xu, K Yan, et al. 2023.  Artificial intelligence for science in quantum, atomistic, and continuum systems. arxiv 2023.  _arXiv preprint arXiv:2307.08423_.
  * Zheng et al. (2024) Chen Zheng, Ke Sun, Hang Wu, Chenguang Xi, and Xun Zhou. 2024.  Balancing enhancement, harmlessness, and general capabilities: Enhancing conversational llms with direct rlhf.  _arXiv preprint arXiv:2403.02513_.
  * Zhong et al. (2023) Ming Zhong, Siru Ouyang, Yizhu Jiao, Priyanka Kargupta, Leo Luo, Yanzhen Shen, Bobby Zhou, Xianrui Zhong, Xuan Liu, Hongxiang Li, et al. 2023.  Reaction miner: An integrated system for chemical reaction extraction from textual data.  In _Proceedings of the 2023 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pages 389–402.

##  Appendix A Language-molecule Translation Instructions

Below is an instruction that describes a task, paired with an input that
provides further context. Write a response that appropriately completes the
request. ### Instruction: You are a researcher. You can come up captions based
on your existing knowledge. Captions are given against the following input.
You should be as detailed as possible. ### Input: Molecule: {source molecule}
In that molecule, could you formulate a caption about? ### Response:{target
caption} Figure 3: Instruction for molecule to language translation, i.e.,
$M\rightarrow L$ Below is an instruction that describes a task, paired with an
input that provides further context. Write a response that appropriately
completes the request. ### Instruction: You are a researcher. You can come up
molecule smile strings based on your existing knowledge. Molecule smile
strings are given against the following input. You should be as detailed as
possible. ### Input: Caption: {source caption} In that caption, could you
generate a molecule smile string? ### Response: {target molecule} Figure 4:
Instruction for language to molecule translation, i.e., $L\rightarrow M$

##  Appendix B Training Effectiveness and Efficiency

Figure 5: Training convergence Figure 6: Training efficiency Figure 7:
Validation loss

##  Appendix C Training Parameters

    args = TrainingArguments(
        output_dir=save_path,
        overwrite_output_dir=True,
        load_best_model_at_end=True,
        num_train_epochs=3,
        per_device_train_batch_size=1
        per_device_eval_batch_size=1
        gradient_accumulation_steps=64
        gradient_checkpointing=False
        optim="adamw_torch_fused",
        learning_rate=5e-5,
        max_grad_norm=0.3,
        warmup_ratio=0.1,
        lr_scheduler_type="cosine",

    )

    (
        load_in_4bit=True,
        bnb_4bit_use_double_quant=True,
        bnb_4bit_quant_type=nf64,
          bnb_4bit_compute_dtype=torch.bfloat16
    )

    (
        lora_alpha=16,
        r = 64,
        lora_dropout=0.1,
        task_type="CAUSAL_LM",
        bias=False,
        target_modules= "all-linear"
    )

Figure 8: Training configuration


# New Cycle-based Formulation, Cost Function, and Heuristics for DC OPF Based
Controlled Islanding

Ilya Tyuryukanov, Marjan Popov, Jorrit A. Bos, Mart A. M. M. van der Meijden,
and Vladimir Terzija I. Tyuryukanov and M. Popov are with Delft University of
Technology, Delft 2628CD, Netherlands. Since October 1, 2021, I. Tyuryukanov
is also with Siemens Energy AG, Freyeslebenstr. 1, Erlangen 91058, Germany
(e-mail<EMAIL_ADDRESS>m.popov@tudelft.nl).J. Bos and M. A. M. M.
van der Meijden are with TenneT TSO B.V., Utrechtseweg 310, Arnhem 6812AR,
Netherlands (e-mail: {Jorrit.Bos, mart.vander.meijden}@tennet.eu). M. A. M. M.
van der Meijden is also with Delft University of Technology, Delft 2628CD,
Netherlands.V. Terzija is with Key Laboratory of Power System Intelligent
Dispatch and Control of Ministry of Education, School of Electrical
Engineering, Shandong University, Jinan 250061, China (email:
<EMAIL_ADDRESS>
###### Abstract

This paper presents a new formulation for intentional controlled islanding
(ICI) of power transmission grids based on mixed-integer linear programming
(MILP) DC optimal power flow (OPF) model. We highlight several deficiencies of
the most well-known formulation for this problem and propose new enhancements
for their improvement. In particular, we propose a new alternative
optimization objective that may be more suitable for ICI than the minimization
of load shedding, a new set of island connectivity constraints, and a new set
of constraints for DC OPF with switching, and a new MILP heuristic to find
initial feasible solutions for ICI. It is shown that the proposed improvements
help to reduce the final optimality gaps as compared to the benchmark model on
several test instances.

###### Index Terms:

Controlled islanding, integer programming, network connectivity, DC OPF,
heuristics

##  I Introduction

Intentional controlled islanding (ICI) is an important system integrity
protection scheme (SIPS) aiming to prevent system collapse due to wide-area
instability by separating the system into a set of non-interacting islands. A
typical situation requiring such a control action is loss of synchronism
between some parts of the system, but it can also be used to limit the spread
of cascading outages, or to isolate unstable parts of the system. The high
relevance of this control action is confirmed by the increasing number of
blackouts around the globe [1].

From the mathematical perspective, computing the splitting boundary to
separate the system into a number of self-sustained islands equates to an OPF
problem with switching. Such problems are notoriously hard to solve, as OPF is
nonlinear, and switching decisions are discrete, which amounts to a mixed-
integer nonlinear program (MINLP). Moreover, ICI must be solved rapidly in
real-time as a corrective control action for instability mitigation. For this
reason, solving ICI as a MINLP is not feasible, and MILP reformulations of the
original MINLP are often seen as an acceptable trade-off between accuracy and
computation time.

Although multiple MILP approximations for ICI and optimal transmission
switching (OTS) exist in literature (e.g., [2, 3]), this paper is focused on
the basic DC OPF ICI model. Our goal is to resolve the fundamental issues
associated with the DC OPF ICI model to devise the efficient computational
enhancements that could also be applied to more advanced MILP-based ICI models
[2, 4, 3].

We first observe that the conventional ICI objective of minimizing the total
post-splitting load shedding [5, 6] may often result in a poor MIP optimality
gap progression if the exact optimal load shedding is close to zero. To
alleviate this issue, we propose to minimize the total load-generation
imbalance in the computed islands while possibly keeping the load shedding as
the secondary optimization objective. Minimizing the total load-generation
imbalance should limit load shedding and additionally limit the value of the
rate of change of frequency (ROCOF) following the system splitting.

Our second contribution aims to improve the modeling of the islands’
connectedness requirement. The previously published papers on optimization-
based ICI either use artificial commodity flows to impose connectivity
constraints on islands [6] or neglect this requirement [7, 2]. Our proposal is
to model connectivity through directed spanning trees, which aims to avoid the
computationally inefficient _big-M_ coefficients that accompany the commodity
flow based formulation. This technique has previously been applied to enforce
radiality in distibution networks [8], but our approach differs from [8].

Our third contribution addresses the big-M coefficients that are normally
present in the DC OPF constraints related to switching decisions. We
reformulate these constraints as Kirchhoff’s voltage law (KVL) equations with
additional slack variables to eliminate the uncertain voltage angle
differences across an open line otherwise modeled as large big-M constants.
Cycle inequalities involving KVL have been previously proposed in [9], but our
work differs from [9] in several aspects. The goal of [9] was to strengthen
the LP relaxation of the conventional DC optimal transmission switching (OTS)
model, while our goal is the elimination of the big-M coefficients. This
motivates a different modeling approach.

Our fourth contribution is a new MILP heuristic to find initial feasible
solutions for DC OPF ICI. For some MIP problems, a feasible solution can be
obtained trivially or in polynomial time, but the OPF ICI problem combining
power flow physics and discrete decisions does not allow for a simple
initialization. Obtaining an early feasible solution is crucial for the
progress of a modern MIP solver, as there are efficient heuristics that are
able to progressively improve the best feasible solution (e.g., the RINS
heuristic). Unfortunately, we could not identify any published MILP heuristics
for OPF ICI, which prompted the design of a simple yet effective method based
on the LP relaxation of DC OPF ICI.

The following sections describe the above contributions in more detail.
Section II-B introduces the used notation and the conventional benchmark DC
OPF ICI model. Section III details our proposed model and cost function.
Section IV introduces the MILP heuristic for OPF ICI and explains how MILP is
solved. Section V describes the test cases and presents the computational
results. Section VI concludes the paper.

##  II Overview of Previous Results

###  II-A Notations

A power network consisting of $n$ nodes and $m$ branches is represented by a
graph with set of nodes $\mathcal{V}$, set of edges $\mathcal{E}$
($|\mathcal{E}|=m$), and set of arcs $\mathcal{A}$ ($|\mathcal{A}|=2m$). Power
demand and power generation at node $i$ prior to splitting are denoted as
$P_{L,i}^{s}$ and $P_{G,i}^{s}$. Load demand and generation shedding to
satisfy the post-splitting conditions are denoted as $P_{L,i}$ and $P_{GS,i}$.
Power flow through branch $(i,j)\in\mathcal{E}$ is defined as $p_{i,j}$, its
value prior to splitting is $p_{i,j}^{s}$, and its limiting magnitude is
$p_{i,j}^{max}$. Admittance of branch $(i,j)$ is denoted as $b_{i,j}$. As the
used model is DC OPF, reactive power is not modeled.

The splitting configuration is modeled by binary variables $x_{i,k}$
($x_{i,k}=1$ implies that node $i$ is assigned to island $k$) and $y_{i,j}$
($y_{i,j}=1$ implies that branch $(i,j)\in\mathcal{E}$ is open). The total
number islands is $K$. An important ICI requirement that promotes transient
stability after splitting is to assign generators from the same coherent group
to one island. Coherent groups are denoted as
$\mathcal{G}_{1},\ldots,\mathcal{G}_{K}$ (i.e., each coherent group is
supposed to form a separate island).

Indices $i$ and $j$ run over network nodes and pairs of nodes (edges or arcs).
Index $k$ runs over islands ($k=1,\ldots,K$). Superscripts $s$, $*$, $LR$
denote pre-splitting values, integral solutions of MILP programs, and linear
relaxation (LR) solutions of MILP programs respectively.

###  II-B Basic MILP Model for DC OPF ICI

Probably the first paper proposing a MILP-based ICI strategy was [5]. Since
then several related models appeared in the literature, of which the model in
[7] is the most compact one. Unlike [5, 10, 11], the model in [7] does not use
auxiliary binary variables that model the membership of edge $(i,j)$ in island
$k$. Instead, the edge status variables $y_{i,j}$ are used directly, which
eliminates $O(mK)$ extra binary variables and constraints. However unlike most
other references, the models in [7, 2] neglect the islands’ connectedness
requirement. To avoid the deficiencies of any single existing model, a
combined DC OPF ICI model is introduced below:

$\displaystyle\text{min}\quad\beta\sum_{i\in\mathcal{V}}P_{LS,i}+\gamma\sum_{i\in\mathcal{V}}P_{GS,i}+\mu\sum_{(i,j)\in\mathcal{E}}p_{i,j}^{s}y_{i,j}$
(1a) $\displaystyle\sum_{k=1}^{K}x_{i,k}=1,~{}\forall i\in\mathcal{V}$ (1b)
$\displaystyle x_{i,k}=1,~{}\forall i\in\mathcal{G}_{k},~{}\forall k$ (1c)
$\displaystyle x_{i,k}-x_{j,k}\leq
y_{i,j},\forall(i,j)\in\mathcal{E},~{}\forall k$ (1d) $\displaystyle
x_{j,k}-x_{i,k}\leq y_{i,j},\forall(i,j)\in\mathcal{E},~{}\forall k$ (1e)
$\displaystyle x_{i,k}+x_{j,k}+y_{i,j}\leq
2,\forall(i,j)\in\mathcal{E},~{}\forall k$ (1f)
$\displaystyle\sum_{(i,j)\in\mathcal{E}}p_{i,j}-\sum_{(j,i)\in\mathcal{E}}p_{j,i}=P_{G,i}^{s}-P_{GS,i}-$
$\displaystyle P_{L,i}^{s}+P_{LS,i},~{}\forall i\in\mathcal{V}$ (1g)
$\displaystyle P_{GS,i}^{min}\leq P_{GS,i}\leq P_{GS,i}^{max},~{}\forall
i\in\mathcal{V}$ (1h) $\displaystyle P_{LS,i}^{min}\leq P_{LS,i}\leq
P_{L,i}^{s},~{}\forall i\in\mathcal{V}$ (1i) $\displaystyle p_{i,j}\leq
p_{i,j}^{max}(1-y_{i,j}),~{}\forall(i,j)\in\mathcal{E}$ (1j) $\displaystyle
p_{i,j}\geq-p_{i,j}^{max}(1-y_{i,j}),~{}\forall(i,j)\in\mathcal{E}$ (1k)
$\displaystyle p_{i,j}-b_{i,j}(\varphi_{i}-\varphi_{j})\leq
M_{i,j}^{\varphi}y_{i,j},~{}\forall(i,j)\in\mathcal{E}$ (1l) $\displaystyle
p_{i,j}-b_{i,j}(\varphi_{i}-\varphi_{j})\geq-
M_{i,j}^{\varphi}y_{i,j},~{}\forall(i,j)\in\mathcal{E}$ (1m)
$\displaystyle\varphi^{min}\leq\varphi_{i}\leq\varphi^{max},~{}\forall
i\in\mathcal{V}$ (1n) $\displaystyle\varphi_{r}=0,~{}r\in\mathcal{R}$ (1o)
$\displaystyle f_{i,j}\leq(n-1)(1-y_{i,j}),~{}\forall(i,j)\in\mathcal{E}$ (1p)
$\displaystyle f_{i,j}\geq-(n-1)(1-y_{i,j}),~{}\forall(i,j)\in\mathcal{E}$
(1q)
$\displaystyle\sum_{(r_{k},j)\in\mathcal{E}}f_{r_{k},j}-\sum_{(j,r_{k})\in\mathcal{E}}f_{j,r_{k}}=\sum_{i\in\mathcal{V}}x_{i,k},~{}$
$\displaystyle r_{k}\in\mathcal{R},~{}k=1,\ldots,K,~{}x_{r_{k},k}=1$ (1r)
$\displaystyle\sum_{(i,j)\in\mathcal{E}}f_{i,j}-\sum_{(j,i)\in\mathcal{E}}f_{j,i}=-1,~{}\forall
i\in\mathcal{V}\setminus\mathcal{R}$ (1s) $\displaystyle
x_{i,k}\in\\{0,1\\},~{}\forall i\in\mathcal{V},~{}k=1,\ldots,K$ (1t)
$\displaystyle y_{i,j}\in\\{0,1\\},~{}\forall(i,j)\in\mathcal{E}$ (1u)

The objective of (1) combines the relevant objectives from [5, 10, 11] with
different weight factors. The factor $\beta$ is associated with minimal load
shedding, which is of great importance for ICI. While minimizing load
shedding, it is important to avoid arbitrary large values of shed or
disconnected generation, which is achieved by introducing the minimization of
generation shedding into (1a) with the weight $\gamma$. Additionally, [10, 11]
use the minimization of total power flow disruption as their cost function. In
our experience, this choice of cost function allows for an easier convergence
to the optimum, but it cannot well enough represent the power balance within
islands, which is the main challenge of ICI. However, a large power flow
disruption is undesirable [12], which explains its presence in (1a) with the
weight $\mu$.

Switching constraints (1b)–(1f) ensure the separation of buses belonging to
different islands from each other through open branches ($y_{i,j}=1$). Here
(1f) ensures that lines connecting buses in the same island are closed, thus
precluding line switching inside islands. Because of (1f), the formulation in
(1) becomes easier to solve, as its feasible region becomes more constrained.
From an operational perspective, requiring both the optimal splitting cutset
and optimal line switching inside of each island may be too complex and
impractical during an emergency condition.

DC OPF constraints (1g)–(1m) describe the physics of DC power flow. Kirchhoff
current law (KCL) is modeled by (1g), which represents the power balance at
each node. Constraints (1j)–(1m) model Ohm’s law under the presence of
switching decisions $y_{i,j}$, which requires the introduction of the node
potential variables $\varphi_{i},~{}\forall i\in\mathcal{V}$. When branch
$(i,j)$ is closed, $p_{i,j}$ is governed by Ohm’s law. If branch $(i,j)$ is
open, $p_{i,j}=0$ and the angle difference $(\varphi_{i}-\varphi_{j})$ can be
arbitrary. Therefore, the big-M constants $M_{i,j}^{\varphi}$ are needed to
allow $(\varphi_{i}-\varphi_{j})$ to take some rather large values.
Additionally, the lower and upper bounds $\varphi^{min}$ and $\varphi^{max}$
in (1n) can also be seen as big-M constants, as their values often cannot be
estimated exactly.

Single commodity flow constraints are used to ensure the connectedness of each
island in the same way as the switching constraints (1b)-(1e) ensure the
separation of islands from one another. For each coherent group
$\mathcal{G}_{k}$, a root node
$r_{k}\in\mathcal{G}_{k},~{}r_{k}\in\mathcal{R}$ is selected, whereby
$\mathcal{R}$ is the union of root nodes of all groups. Constraint (1r)
requires that the number of units of an artificial commodity produced at
$r_{k}$ equals the total number of nodes in island $k$ corresponding to
$\mathcal{G}_{k}$, while constraints (1s) demand one unit of commodity to be
consumed at each non-root node. The satisfaction of flow balance (1r)–(1s)
requires each island to be connected. However, the upper bounds of commodity
flows for each branch are not known, which prompts the usage of another big-M
constant $n-1$ in (1p)–(1q). This big-M constant could be somewhat lowered
through exact calculations, but $n-1$ is shown here to simplify the
representation. Unlike [5, 10], (1) does not require a separate set of
commodity flow variables for each island, as it appears that a single set of
$m$ commodity flows $f_{ij}$ can be used for all islands.

##  III Proposed MILP Model

###  III-A DC OPF ICI using Cycle-based KVL Constraints

A major issue with (1) is the presence of big-M constants $M_{i,j}^{\varphi}$
in (1l)–(1m). While there are methods to bound these coefficients when the
solution is required to contain a single connected component (e.g., [13]), no
bound strengthening method is known for the disconnected networks resulting
from ICI. Thus, a possibly unrealistically large value needs to be assumed for
$M_{i,j}^{\varphi}$ in order to preserve all feasible solutions.

An alternative to Ohm’s law form of power flow equations in (1l)–(1m) is to
apply KVL to the entire network, while considering the switching decisions.
The KVL equations are written with respect to a set of linearly independent
cycles in the network, which is known as _cycle basis_. For the reasons that
become clear below, it is desirable to compute the cycle basis as quickly as
possible, which is achieved by using as cycle basis _fundamental cycles_
induced by the network’s minimal spanning tree (MST) [14]. Given cycle basis
$\mathcal{CB}$, the KVL alternatives of (1l)–(1m) can be formulated as
follows:

$\displaystyle\sum_{(i,j)\in
C}\frac{p_{i,j}}{b_{i,j}}\text{sgn}(i,j)\leq\frac{1}{2}M_{C}\sum_{(i,j)\in
C}y_{i,j},~{}\forall C\in\mathcal{CB}$ (2a) $\displaystyle\sum_{(i,j)\in
C}\frac{p_{i,j}}{b_{i,j}}\text{sgn}(i,j)\geq-\frac{1}{2}M_{C}\sum_{(i,j)\in
C}y_{i,j},~{}\forall C\in\mathcal{CB}$ (2b)

where $C\in\mathcal{CB}$ is a cycle formed by undirected edges
$(i,j)\in\mathcal{E}$, and function $\text{sgn}(i,j)$ takes the plus sign if
edge $(i,j)$ is aligned with the clockwise direction of $C$ and the minus sign
otherwise.

As (1f) precludes line switching inside islands, each cycle is either
connected or has at least two open lines. Therefore, $M_{C}$ is the sum of
$|\frac{p_{i,j}^{max}}{b_{i,j}}|,(i,j)\in C$ minus the two smallest values of
$|\frac{p_{i,j}^{max}}{b_{i,j}}|$ along $C$. Unlike $M_{i,j}^{\varphi}$,
$M_{C}$ can be computed exactly, but it may still be quite large for long
cycles. Thus, to additionally strengthen the model, we are adding some
additional KVL constraints (2) associated with cycles shorter than a certain
length. All such cycles can be found in polynomial time by applying to each
node techniques based on recursive node neighborhood traversal. Short cycles
(e.g., up to length 7) can be found quickly before the start of the MIP
solver, but enumerating longer cycles becomes increasingly computationally
inefficient.

Obviously, any cycle basis of the disconnected network contains less cycles
than the cycle basis of the original network. As only a small fraction of
edges is opened for ICI, the majority of fundamental cycles of the
disconnected network coincide with the initial fundamental cycles. However, in
some cases the original fundamental cycles do not fully describe the cycle
basis of the disconnected network. Thus, for each integer solution obtained
during the solution process, cycles violating (2) need to be identified, and
the corresponding inequalities (2) added to the model by using MIP solver
callbacks. A possible computationally efficient way to achieve this is as
follows:

  1. 1.

Compute the MST of the original network and the associated fundamental cycle
basis $\mathcal{CB}_{0}$.

  2. 2.

Assign close-to-zero weights to the network edges belonging to the MST.

  3. 3.

Once an integer solution is found, compute its minimal spanning forest and the
associated fundamental cycle basis $\mathcal{CB}_{i}$ ($i$ is the solution
number). The previously assigned small edge weights should promote the
alignment of the newly computed cycle basis with the original one.

  4. 4.

Check the solution for conformity with (2) (i.e., $\sum_{(i,j)\in
C}\frac{p^{*}_{i,j}}{b_{i,j}}\text{sgn}(i,j)=0,~{}\forall
C\in\mathcal{CB}_{i}$).

  5. 5.

If any cycles in $\mathcal{CB}_{i}$ violate (2), use them to add new KVL
constraints (2) to the model.

Besides strengthening the model, the abovementioned KVL constraints based on
extra short cycles also improve the satisfaction of (2) for the intermediate
integer feasible solutions.

###  III-B Island Connectivity Constraints

The single commodity flow constraints that are most commonly used to enforce
islands’ contiguity [5, 10, 6] have the drawback of introducing at least $m$
auxiliary variables and large big-M constants that may lead to loose LP
relaxations. In [11], two alternative methods based on multicommodity flows
and network cutsets are proposed. However, these alternatives may often
underperform the simple single commodity flow approach. At least the
multicommodity flow method has largely theoretical importance due to the very
large number of auxiliary variables that it requires [15]. The cutset-based
method requires solving many max-flow min-cut problems to identify cutsets
that violate connectivity, which has the time complexity of $O(mn^{2})$ for a
single max-flow min-cut run.

In this section, a different approach to islands’ connectivity is proposed
that does not involve big-M coefficients and allows to quickly find strong
inequalities that ensure network connectivity. This approach is related to the
studies on radiality constraints in distribution networks (e.g., [16, 8]), but
has some application-specific features and computational improvements. One of
conditions for network connectedness is that there exists a spanning tree with
$n-1$ edges that spans all of its nodes. In general, for a network with $k$
connected components there exists $k$ spanning trees (one for each component)
forming a _spanning forest_ that includes $n-k$ edges. This observation
motivates the following connectivity formulation based on rooted directed
spanning trees:

$\displaystyle\sum_{(i,j)\in\mathcal{A}}z_{i,j}=n-k$ (3a) $\displaystyle
z_{j,r}=0,~{}\forall(j,r)\in\mathcal{A},~{}\forall r\in\mathcal{R}$ (3b)
$\displaystyle\sum_{(r,j)\in\mathcal{A}}z_{r,j}\geq 1,~{}\forall
r\in\mathcal{R}$ (3c)
$\displaystyle\sum_{(j,i)\in\mathcal{A}}z_{j,i}=1,~{}\forall
i\in\mathcal{V}\setminus\mathcal{R}$ (3d) $\displaystyle z_{i,j}+z_{j,i}\leq
1-y_{i,j},~{}\forall(i,j)\in\mathcal{E}_{C}$ (3e) $\displaystyle
z_{i,j}+z_{j,i}=1-y_{i,j},~{}\forall(i,j)\in\mathcal{E}\setminus\mathcal{E}_{C}$
(3f) $\displaystyle\sum_{(i,j)\in C}z_{i,j}\leq|C|-1,~{}\forall
C\in\mathcal{D}$ (3g) $\displaystyle
z_{i,j}\in\\{0,1\\},~{}\forall(i,j)\in\mathcal{A}$ (3h)

where $z_{i,j}$ and $z_{j,i}$ represent the status of arcs $(i,j)$ and
$(j,i)\in\mathcal{A}$ ($z_{i,j}=1$ iff $(i,j)$ is enabled).

The equalities in (3b)–(3c) define the root nodes of $K$ coherent generator
groups as root nodes of $K$ trees each spanning one island. According to
(3b)–(3c), each root node should have no incoming arcs and at least one
outgoing arc. Conversely, each non-root node should have one incoming arc
according to (3d). Constraints (3e)–(3f) link $z$ and $y$ variables. For the
edges that participate in network cycles (i.e., in the cycle basis described
in Section III-A), (3e) is valid. For the rest of the edges, the stronger
version (3f) can be used instead of (3e). The number of enabled arcs is
constrained by (3a).

Combined constraints (3a)–(3f) can define $k$ connected components spanned by
$k$ directed trees, or more than $k$ connected components. In the latter case,
the additional components will form directed cycles to satisfy (3a)–(3f).
Therefore, cycle breaking constraints (3g) are needed to only allow spanning
forests as solutions, where $\mathcal{D}$ is the set of all directed cycles in
the network and $C$ is a sequence of arcs in $\mathcal{A}$ forming a directed
cycle. As there are exponentially many cycles in the network, all constraints
in (3g) cannot be simultaneously added to the model. Instead, only the
violated inequalities (3g) are added in MIP solver callbacks that are
available in most of modern MIP solvers (see Section III-E).

###  III-C Extra Strengthening Cycle Inequalities

To further strengthen the cycle-based ICI formulation, the following _triangle
inequalities_ (e.g., see [17]) are added for all known cycles of the minimal
length:

$y_{i,j}-\sum_{\begin{subarray}{c}(p,q)\in C\\\
(p,q)\neq(i,j)\end{subarray}}y_{p,q}\leq 0,~{}\forall
C\in\mathcal{CB},|C|=3,\forall(i,j)\in C$ (4)

###  III-D Minimization of Load-Generation Imbalance

Minimal load shedding is often chosen as the cost function of ICI, as it is
highly important for utilities to minimize the loss of load. However, shedding
an exact amount of load at exact bus at exact moment in time to realize the
outcome of (1) may be sometimes quite challenging with the available
equipment. Moreover, the optimal load shedding can be zero, in which case it
will be impossible to provide a non-zero lower bound, and the integrality gap
of the MIP solver will be close to 100 % during the whole optimization.

Among practitioners, the notion of _balanced islands_ with good load-
generation balance is no less popular than the idea of minimal load shedding.
Islands with good power balance would naturally limit load shedding without
the need of picking specific buses, while also limiting the ROCOF value after
splitting. At the same time, the total power imbalance of all islands is
always non-zero and can be much easier bounded from below. Based on (1), the
minimization of load generation imbalance can be formulated as follows:

$\displaystyle\text{min}\quad\sum_{k=1}^{K}P_{\Delta,k}$ (5a)
$\displaystyle\quad
P_{\Delta,k}\geq\sum_{i\in\mathcal{V}}(P_{G,i}^{s}-P_{L,i}^{s})x_{i,k},~{}\forall
k$ (5b) $\displaystyle
P_{\Delta,k}\geq\sum_{i\in\mathcal{V}}(P_{L,i}^{s}-P_{G,i}^{s})x_{i,k},~{}\forall
k$ (5c) $\displaystyle\eqref{eqn:dc_xik}-\eqref{eqn:XYbin}$

where $P_{\Delta,k}$ represents the power imbalance of island $k$.

The inequalities in (5b)–(5c) are constraining the absolute value
$|\sum_{i\in\mathcal{V}}(P_{G,i}^{s}-P_{L,i}^{s})x_{i,k}|$ from below for each
$k$, while the minimization requirement (5a) constrains it from above.

###  III-E Overall Solution Process

The ideas from Sections III-B–III-A can be combined together into the
following novel DC OPF ICI formulation:

$\displaystyle\text{min}~{}\alpha\sum_{k=1}^{K}P_{\Delta,k}+\beta\sum_{i\in\mathcal{V}}P_{LS,i}+\gamma\sum_{i\in\mathcal{V}}P_{GS,i}+\mu\sum_{(i,j)\in\mathcal{E}}p_{i,j}^{s}y_{i,j}$
$\displaystyle\eqref{eqn:dc_xik}-\eqref{eqn:dc_Pij_Lo}$
$\displaystyle\eqref{eqn:KVL_Y_Hi}-\eqref{eqn:KVL_Y_Lo}$ (6)
$\displaystyle\eqref{eqn:sptree}-\eqref{eqn:z01},~{}\eqref{eqn:triang},~{}\eqref{eqn:genimb}-\eqref{eqn:lodimb}$

The above formulation requires MILP callbacks for cycle breaking constraints
(3g) to be solved efficiently. In our implementation, the violated cycle
breaking constraints are easily separated using Algorithm 1 for each new
integer solution.

1:Input: $z_{i,j}^{*}$, $\mathcal{V}$, $K$

2:Output: Set of violating directed cycles $\mathcal{C}$

3:$G(\mathcal{V},\mathcal{A}_{G})\leftarrow\text{directed graph induced
by~{}}z_{i,j}^{*}$

4:$F(\mathcal{V},\mathcal{A}_{F})\leftarrow\mathtt{spanning\\_forest}(G)$

5:$\mathcal{CC}\leftarrow\mathtt{connected\\_components}(F)$

6:if $|\mathcal{CC}|=K$ then return $\emptyset$

7:end if

8:$\mathcal{A}_{C}\leftarrow\mathcal{A}_{G}\setminus\mathcal{A}_{F}$

9:for $a_{i,j}\in\mathcal{A}_{C}$ do

10:
$P\leftarrow\mathtt{shortest\\_path}(\mathcal{F},\text{source=}j,\text{target=}i)$

11: $C\leftarrow P\cup a_{i,j}$

12: $\mathcal{C}\leftarrow\mathcal{C}\cup C$

13:end for

14:return $\mathcal{C}$

Algorithm 1 Generation of cycle breaking constraints

Given an integer solution $z_{i,j}^{*}$, it can be considered as directed
graph $G$. For this graph, the spanning forest $F$ can be computed by using
the standard spanning tree algorithms with loglinear time complexity. If $F$
has $k$ connected components, the solution is accepted, and the best integer
solution gets updated. Otherwise, Algorithm 1 finds arcs that are in $G$, but
not in $F$. Each such arc corresponds to a directed cycle, which allows to
retrieve the full set of cycles $\mathcal{C}$ that violate (3g) in loglinear
time.

In addition to cycle breaking constraints (3g), we are adding directed cutset
constraints between the root node of each group and the remaining generators
of the same group using the same $z_{i,j}$ variables. These constraints are
convenient for enforcing connectivity of a defined set of _terminal nodes_ to
the root node. A detailed description of these constraints including
improvements like _back cuts_ and _creep flow_ can be found in [15] and
several other references.

The cycle breaking constraints in (3g) also combine well with the cycle-based
DC OPF ICI formulation in (2). To apply (2), a cycle basis and optionally some
additional network cycles need to be found. From these precomputed cycles, the
shorter ones (e.g., with the length less than 6–8 edges) can be used to
strengthen the initial model formulation with some constraints of type (3g).

##  IV MILP Heuristics

When solving (6), no feasible solution could be found within the optimization
time limit for several larger instances. In practice, it is highly beneficial
to obtain an initial feasible solution as early as possible because this would
increase the available time to improve the upper bound of (6) using various
efficient techniques available in the modern MIP solvers. As the practice has
shown that multiple MILP heuristics available in Gurobi [18] often cannot
retrieve a feasible solution, a new problem-specific heuristic had to be
proposed to improve the solvability of (1) and (6). This heuristic is outlined
in Algorithm 2.

1:Input: $x_{i,k}^{LP}$, $\mathcal{V}$, $\mathcal{E}$, $\mathcal{R}$, K

2:Output: $x_{i,k}^{*}$, $y_{i,j}^{*}$, $z_{i,j}^{*}$

3:$\mathcal{E}^{on}\leftarrow\emptyset$

4:for $k=1,\ldots,K$ do

5: for $(i,j)\in\mathcal{E}$ do

6: if $x_{i,k}^{LP}>0.9\wedge x_{j,k}^{LP}>0.9$ then
$\mathcal{E}^{on}\leftarrow\mathcal{E}^{on}\cup(i,j)$

7: end if

8: end for

9:end for

10:$\mathcal{CC}\leftarrow\mathtt{connected\\_components}(G(\mathcal{V},\mathcal{E}^{on}))$

11:$\mathcal{P}\leftarrow\emptyset$ // Set of partial islands

12:for $CC\in\mathcal{CC}$ do

13: $r\leftarrow\mathcal{R}\cap CC$

14: if $|r|=0$ then continue

15: else if $|r|\geq 2$ then return $\emptyset$

16: else$~{}\mathcal{P}\leftarrow\mathcal{P}\cup CC$

17: end if

18:end for

19:if $\frac{\sum_{P\in\mathcal{P}}|P|}{n}\geq 0.8$ then

20: Fix $x_{i,k}$, $y_{i,j}$, $z_{i,j}$ according to $\mathcal{P}$ and solve
the reduced MILP (e.g., (1)) with more than 80% fixed buses.

21:end if

22:return The solution of reduced MILP $x_{i,k}^{*}$, $y_{i,j}^{*}$,
$z_{i,j}^{*}$

Algorithm 2 MILP heuristic for initialization of DC OPF ICI

The LR of $x_{i,k}$ is taken as input and analyzed in lines 4–9 to identify
the edges that are considered to be closed by the current LR. Then these edges
are used to build an undirected graph $G$ in line 10, and the connected
components of that graph are identified. Next, each connected component is
analyzed for the presence of a root node $r\in\mathcal{R}$ in it. If a root
node is present is the component, it is saved to the set of _partial islands_
$\mathcal{P}$. Finally, the $x_{i,k}$ variables of the nodes belonging to
contiguous partial islands are fixed, and the values of the remaining node
variables are obtained by solving the initial problem as MILP with the
majority of discrete variables being fixed. Because of a high degree of
reduction of discrete variables, the reduced MILP is usually solved in a
fraction of a second. However, it is important to find connected partial
islands, as just rounding the LR will very often result in infeasibility.

The heuristic in Algorithm 2 is only run until a feasible solution becomes
available, as modern MILP solvers are able to apply more efficient heuristics
to consistently improve the existing solution.

##  V Computational Results

###  V-A Test Setup

This section111The code implementation can be found at
https://github.com/ityuryukanov/power_network_clustering under
.papers/PSCC_2022. presents a comparison of the proposed formulation in (6) to
the known benchmark model in (1). The results are produced for a number of
test networks from the MATPOWER toolbox [19]. As publicly available power
system test cases do not include coherent generator groupings that are
required in (1c), those were obtained with the generator coherency algorithm
from [20] for every test case. For the MATPOWER test cases lacking the dynamic
generator data it was fitted from [21].

Three test cases from MATPOWER were selected to evaluate the proposed
formulation: _case89pegase_ , _case1354pegase_ , and _case1888rte_. The 89 bus
test case was selected because it posed some difficulties for the MILP solver
despite its moderate size, while the 1354 and 1888 bus test cases are examples
of large-scale networks.

All computations were performed using Gurobi v. 9.1. [18] on a PC with an
Intel® Core™ i7 2.20 GHz CPU with 6 cores and 16 Gb of RAM running on Windows
10. The reported solution times include the presolve and MILP solver time, but
exclude the model preparation time, which is polynomial. The time limit for
the test cases with less than 500 buses has been set to 480 seconds and to 720
seconds for the test cases with more than 500 buses. In all the experiments,
the default Gurobi solver settings were used. The maximal time to be spent in
the MILP heuristic from Section IV is set to 3 % of the time limit. For
simplicity, the auxiliary MILP model in Algorithm 2 is based on (1) for all
test cases. The solution was considered optimal if its MIP optimality gap was
less than 1%.

TABLE I: Load Shedding Minimization using (1) $n$ | $K$ |  $UB$, | $g,~{}\%$ |  $P_{LS}$, |  $P_{\Delta}$, |  $P_{GS}$, |  $p_{ij}^{\Sigma}$,
---|---|---|---|---|---|---|---
|  | p.u. |  $T$, s | p.u. | p.u. | p.u. | p.u.
89 | 2 | 6.4089 | 47 % | 5.02 | 11.4 | 6.4 | 13.3
89 | 3 | 8.3011 | 46 % | 5.02 | 11.4 | 6.4 | 32.2
89 | 4 | 13.844 | 64 % | 10.2 | 21.9 | 11.6 | 34.9
89 | 5 | 13.876 | 70 % | 10.2 | 21.9 | 11.6 | 35.2
1354 | 2 | 2.3878 | 0.5 s | -0.00 | 16.6 | 16.6 | 22.2
1354 | 3 | 4.6645 | 100 s | -0.00 | 16.6 | 16.6 | 45
1354 | 4 | 3.3921 | 1 s | -0.00 | 16.6 | 16.6 | 32.3
1354 | 5 | 8.3761 | 20 % | -0.00 | 16.6 | 16.6 | 82.1
1354 | 6 | 12 | 12 % | 0.03 | 16.7 | 16.7 | 118
1354 | 7 | 13.631 | 14 % | 0.01 | 16.7 | 16.6 | 135
1354 | 8 | 21.148 | 33 % | 1.58 | 19.8 | 18.2 | 194
1888 | 2 | 0.4949 | 0.8 s | -0.00 | 9.81 | 9.81 | 3.97
1888 | 3 | 11.38 | 51 % | 0.08 | 9.81 | 9.88 | 112
1888 | 4 | 4.5363 | 11 % | -0.00 | 9.81 | 9.81 | 44.4
1888 | 5 | 7.0738 | 14 % | -0.00 | 9.81 | 9.81 | 69.8
1888 | 6 | 9.2755 | 10 % | 0.02 | 9.85 | 9.83 | 91.5
1888 | 7 | 16.976 | 40 % | 0.17 | 9.94 | 9.97 | 167
1888 | 8 | 12.457 | 17 % | -0.00 | 9.81 | 9.81 | 124

###  V-B Minimization of Load Shedding

The first set of test results is related to minimizing the known objective in
(1a). The objective weights are set as follows: $\alpha=0$, $\beta=1$,
$\gamma=0.01$, $\mu=0.1$. The relatively high value of $\mu$ is selected
mainly to ensure non-zero lower bounds in all cases and thus a more
informative progression of the MIP optimality gap during the optimization. The
big-M coefficients $M_{i,j}^{\varphi}$ in (6) are set to $2\pi$,
$\varphi^{min}$ is set to $-\pi$, $\varphi^{max}$ is set to $\pi$ , and the
power flow limitation $p_{ij}^{max}$ is assumed to be equal to
$b_{i,j}\pi/4,~{}\forall(i.j)\in\mathcal{E}$.

First, the results for the baseline formulation in (1) are presented in Table
I. The table header can be read as follows: $UB$ is the best MIP objective
value (i.e., the upper bound), $g$ is the final optimality gap, $T$ is the
solution time, $P_{LS}$, $P_{\Delta}$, $P_{GS}$, $p_{ij}^{\Sigma}$ are the
final values of load shedding, total power imbalance, generator shedding, and
power flow cut respectively. In Table I, $g$ and $T$ are included
interchangeably: if $g$ is less than the optimality tolerance then $T$ is
included, otherwise $g$ is shown, as $T$ equals to the time limit.

Next, in the formulation in (1) the flow-based connectivity constraints
(1p)–(1s) are replaced with (3) and the MILP heuristic is enabled. The results
obtained after this modification are summarized in Table II. Finally, the
results of our complete formulation in (6) with the MILP heuristic enabled for
the given cost function are given in Table III.

As it can be seen, the results in Tables II–III in terms of the MIP optimality
gap and solution time usually exceed those in Table I, which demonstrates the
efficacy of the proposed improvements.

TABLE II: Load Shedding Minimization using (1) and (3) $n$ | $K$ |  $UB$, | $g,~{}\%$ |  $P_{LS}$, |  $P_{\Delta}$, |  $P_{GS}$, |  $p_{ij}^{\Sigma}$,
---|---|---|---|---|---|---|---
|  | p.u. |  $T$, s | p.u. | p.u. | p.u. | p.u.
89 | 2 | 6.4089 | 49 % | 5.02 | 11.4 | 6.4 | 13.3
89 | 3 | 8.3011 | 46 % | 5.02 | 11.4 | 6.4 | 32.2
89 | 4 | 13.844 | 65 % | 10.2 | 21.9 | 11.6 | 34.9
89 | 5 | 13.876 | 68 % | 10.2 | 21.9 | 11.6 | 35.2
1354 | 2 | 2.3878 | 0.5 s | -0.00 | 16.6 | 16.6 | 22.2
1354 | 3 | 4.6645 | 80 s | -0.00 | 16.6 | 16.6 | 45
1354 | 4 | 3.3921 | 1.1 s | -0.00 | 16.6 | 16.6 | 32.3
1354 | 5 | 8.1682 | 16 % | -0.00 | 16.6 | 16.6 | 80
1354 | 6 | 12 | 12 % | 0.03 | 16.7 | 16.7 | 118
1354 | 7 | 13.453 | 13 % | 0.03 | 16.7 | 16.7 | 133
1354 | 8 | 20.278 | 30 % | 0.29 | 17.2 | 16.9 | 198
1888 | 2 | 0.4949 | 1 s | -0.00 | 9.81 | 9.81 | 3.97
1888 | 3 | 9.0864 | 39 % | 0.01 | 9.82 | 9.81 | 89.8
1888 | 4 | 4.5119 | 10 % | -0.00 | 9.81 | 9.81 | 44.1
1888 | 5 | 7.0672 | 14 % | -0.00 | 9.81 | 9.81 | 69.7
1888 | 6 | 9.3695 | 12 % | -0.00 | 9.81 | 9.81 | 92.7
1888 | 7 | 14.329 | 28 % | -0.00 | 9.81 | 9.81 | 142
1888 | 8 | 12.158 | 15 % | 0.00 | 9.81 | 9.81 | 121
TABLE III: Load Shedding Minimization using (6) $n$ | $K$ |  $UB$, | $g,~{}\%$ |  $P_{LS}$, |  $P_{\Delta}$, |  $P_{GS}$, |  $p_{ij}^{\Sigma}$,
---|---|---|---|---|---|---|---
|  | p.u. |  $T$, s | p.u. | p.u. | p.u. | p.u.
89 | 2 | 6.4089 | 57 % | 5.02 | 11.4 | 6.4 | 13.3
89 | 3 | 8.3011 | 57 % | 5.02 | 11.4 | 6.4 | 32.2
89 | 4 | 11.254 | 64 % | 5.9 | 13.2 | 7.28 | 52.9
89 | 5 | 11.096 | 62 % | 5.9 | 13.2 | 7.28 | 51.3
1354 | 2 | 2.3878 | 0.5 s | -0.00 | 16.6 | 16.6 | 22.2
1354 | 3 | 4.6645 | 55 s | -0.00 | 16.6 | 16.6 | 45
1354 | 4 | 3.3921 | 1.1 s | -0.00 | 16.6 | 16.6 | 32.3
1354 | 5 | 8.1682 | 17 % | -0.00 | 16.6 | 16.6 | 80
1354 | 6 | 12 | 11 % | 0.03 | 16.7 | 16.7 | 118
1354 | 7 | 13.441 | 13 % | -0.00 | 16.6 | 16.6 | 133
1354 | 8 | 18.420 | 21 % | 0.24 | 17.1 | 16.9 | 180
1888 | 2 | 0.4949 | 0.6 s | -0.00 | 9.81 | 9.81 | 3.97
1888 | 3 | 8.5821 | 33 % | -0.00 | 9.81 | 9.81 | 84.9
1888 | 4 | 4.5116 | 590 s | -0.00 | 9.81 | 9.81 | 44.1
1888 | 5 | 7.0534 | 9 % | -0.00 | 9.81 | 9.81 | 69.6
1888 | 6 | 9.275 | 11 % | 0.02 | 9.85 | 9.83 | 91.5
1888 | 7 | 14.163 | 26 % | -0.00 | 9.81 | 9.81 | 141
1888 | 8 | 12.16 | 13 % | -0.00 | 9.81 | 9.81 | 121

###  V-C Minimization of Total Power Imbalance

To complement Section V-B, the same test cases are recomputed with the main
objective of power imbalance minimization. The objective weights are set as
follows: $\alpha=1$, $\beta=0.01$, $\gamma=0.01$, $\mu=0.01$. Now the highest
priority is given to generation-load imbalance minimization, with load
shedding being included with a small weight. As the initial power imbalance is
nearly always substantially higher than zero after the network is split, non-
zero lower bounds can be produced without assigning an increased value to
$\mu$.

The results comparison follows the same route as in Section V-B. The results
of minimizing the new objective with the formulation in 1 are listed in Table
IV. The outcome of the exchange of (1p)–(1s) with (3) is illustrated in Table
V. Finally, the results of our formulation in (6) are shown in Table VI.

As it can be seen from Table IV, with the change of the objective function
more cases fail to compute at least a feasible solution with the basic
formulation (1). Thus, the role of feasibility MILP heuristics becomes larger.
In fact, the proposed MILP heuristic is enabled for the results in Tables
V–VI, and feasibility is achieved for all the test cases. The performance is
also improved with the proposed enhancements, which can be seen by comparing
the 4th columns of Tables IV–VI.

An additional important observation can be made by comparing the upper bounds
of _case89pegase_ across Tables IV–VI. Although the optimal solution could be
reached in all cases, the optimal value is better with our complete
formulation (6). This is explained by the absence of big-M constraints in (6),
the values for which must be selected. Thus, (6) can be certain to achieve the
true optimal solution. In the present case study, it is possible to achieve
the same optimum using (1) if the angular big-M constraints are substantially
increased (e.g., by a factor of 10).

TABLE IV: Islands’ Power Imbalance Minimization using (1) $n$ | $K$ |  $UB$, | $g,~{}\%$ |  $P_{LS}$, |  $P_{\Delta}$, |  $P_{GS}$, |  $p_{ij}^{\Sigma}$,
---|---|---|---|---|---|---|---
|  | p.u. |  $T$, s | p.u. | p.u. | p.u. | p.u.
89 | 2 | 7.8621 | 5 s | 3.09 | 6.24 | 4.47 | 55.5
89 | 3 | 8.0531 | 33 s | 3.09 | 6.24 | 4.47 | 74.6
89 | 4 | 14.807 | 13 s | 8.77 | 11.1 | 10.1 | 66.6
89 | 5 | 14.809 | 16 s | 8.77 | 11.1 | 10.1 | 66.9
1354 | 2 | 17.023 | 0.5 s | -0.00 | 16.6 | 16.6 | 22.2
1354 | 3 | 17.272 | 11 s | -0.00 | 16.6 | 16.6 | 47.1
1354 | 4 | 17.124 | 1.7 s | -0.00 | 16.6 | 16.6 | 32.3
1354 | 5 | 17.684 | 580 s | -0.00 | 16.6 | 16.6 | 88.3
1354 | 6 | – | – | – | – | – | –
1354 | 7 | – | – | – | – | – | –
1354 | 8 | – | – | – | – | – | –
1888 | 2 | 9.9447 | 1.2 s | -0.00 | 9.81 | 9.81 | 3.97
1888 | 3 | – | – | – | – | – | –
1888 | 4 | 10.427 | 23 s | -0.00 | 9.81 | 9.81 | 52.2
1888 | 5 | 10.603 | 94 s | -0.00 | 9.81 | 9.81 | 69.8
1888 | 6 | 10.896 | 120 s | -0.00 | 9.81 | 9.81 | 99.1
1888 | 7 | – | – | – | – | – | –
1888 | 8 | – | – | – | – | – | –
222This work was financially supported by the Dutch Scientific Council NWO in
collaboration with TSO TenneT, DSOs Alliander, Stedin, Enduris, VSL and
General Electric in the framework of the Energy System Integration & Big Data
program under the project ”Resilient Synchromeasurement-based Grid Protection
Platform, no. 647.003.004”. The authors also greatly thank the anonymous
reviewers for their highly useful suggestions.

##  VI Conclusions

This paper has proposed several computational enhancements to the DC OPF ICI
problem. The main idea was to eliminate the big-M constants that are present
in the most existing formulations in order to tighten the linear relaxation of
DC OPF ICI. The big-M constants present in island connectivity constraints
have been removed by replacing the popular single commodity flow based
connectivity model with the new model based on directed spanning forests,
which does not require any large coefficients. The big-M constants associated
with ICI switching decisions and DC OPF in general have been excluded by
replacing the DC OPF switching constraints based on Ohm’s law with the new
ones based on KVL.

Besides of handling of all the three types of big-M constants in the standard
DC OPF ICI formulation, it has also been observed that the existing ICI
objective functions may not fully correspond to the practical requirements
associated with system splitting. To ameliorate this, the new objective
function minimizing the load-generation imbalance after splitting has been
introduced. With this objective, situations when the MIP solver cannot find an
initial feasible solution have become more common. To fix this, a new MILP
heuristic for DC OPF ICI has been proposed as well. In general, the solution
process was noticeably influenced by the type of objective function. Thus, it
could be possible to choose ”convenient” objective functions to quickly find
feasible initial solutions.

Finally, the proposed enhancements have been tested against the compact and
efficient baseline DC ICI OPF model from the literature, and it has been found
that the novel ICI model produces better results in most of cases. In the
future, it is planned to apply the findings of this paper to more complex and
accurate ICI models and possibly to other transmission switching problems.

TABLE V: Power Imbalance Minimization using (1) and (3) $n$ | $K$ |  $UB$, | $g,~{}\%$ |  $P_{LS}$, |  $P_{\Delta}$, |  $P_{GS}$, |  $p_{ij}^{\Sigma}$,
---|---|---|---|---|---|---|---
|  | p.u. |  $T$, s | p.u. | p.u. | p.u. | p.u.
89 | 2 | 7.8621 | 13 s | 3.09 | 6.24 | 4.47 | 55.5
89 | 3 | 8.0531 | 85 s | 3.09 | 6.24 | 4.47 | 74.6
89 | 4 | 14.807 | 13 s | 8.77 | 11.1 | 10.1 | 66.6
89 | 5 | 14.809 | 15 s | 8.77 | 11.1 | 10.1 | 66.9
1354 | 2 | 17.023 | 0.6 s | -0.00 | 16.6 | 16.6 | 22.2
1354 | 3 | 17.32 | 19 s | -0.00 | 16.6 | 16.6 | 51.9
1354 | 4 | 17.123 | 1.5 s | -0.00 | 16.6 | 16.6 | 32.3
1354 | 5 | 17.72 | 90 s | 0.07 | 16.6 | 16.7 | 89.6
1354 | 6 | 18.001 | 220 s | -0.00 | 16.6 | 16.6 | 120
1354 | 7 | 18.162 | 410 s | -0.00 | 16.6 | 16.6 | 136
1354 | 8 | 18.628 | 1.1 % | 0.00 | 16.6 | 16.6 | 183
1888 | 2 | 9.9447 | 1 s | -0.00 | 9.81 | 9.81 | 3.97
1888 | 3 | 11.135 | 3.5 % | -0.00 | 9.81 | 9.81 | 123
1888 | 4 | 10.374 | 36 s | -0.00 | 9.81 | 9.81 | 46.9
1888 | 5 | 10.688 | 280 s | -0.00 | 9.81 | 9.81 | 78.3
1888 | 6 | 10.897 | 280 s | 0.02 | 9.85 | 9.83 | 93.8
1888 | 7 | 11.301 | 170 s | -0.00 | 9.81 | 9.81 | 140
1888 | 8 | 11.211 | 1.1 % | -0.00 | 9.81 | 9.81 | 131

## References

  * [1] H. H. Alhelou, M. E. Hamedani-Golshan, T. C. Njenda, and P. Siano, “A survey on power system blackout and cascading events: Research motivations and challenges,” _Energies_ , vol. 12, no. 4, p. 682, 2019.
  * [2] P. A. Trodden, W. A. Bukhsh, A. Grothey _et al._ , “Optimization-based islanding of power networks using piecewise linear ac power flow,” _IEEE Trans. Power Syst._ , vol. 29, no. 3, pp. 1212–1220, 2014.
  * [3] W. E. Brown and E. Moreno-Centeno, “Transmission-line switching for load shed prevention via an accelerated linear programming approximation of ac power flows,” _IEEE Trans. Power Syst._ , vol. 35, no. 4, pp. 2575 – 2585, 2020\.
  * [4] C. Coffrin and P. Hentenryck, “Transmission system restoration with co-optimization of repairs, load pickups, and generation dispatch,” _Int. J. Electr. Power Energy Syst._ , vol. 72, pp. 144 – 154, 2015.
  * [5] N. Fan, D. Izraelevitz, F. Pan, P. M. Pardalos, and J. Wang, “A mixed integer programming approach for optimal power grid intentional islanding,” _Energy Syst._ , vol. 3, pp. 77–93, 2012.
  * [6] F. Teymouri, T. Amraee, H. Saberi, and F. Capitanescu, “Toward controlled islanding for enhancing power grid resilience considering frequency stability constraints,” _IEEE Trans. Smart Grid_ , vol. 10, no. 2, pp. 1735–1746, 2017\.
  * [7] P. A. Trodden, W. A. Bukhsh, A. Grothey, and K. I. M. McKinnon, “Milp formulation for controlled islanding of power networks,” _Int. J. Electr. Power Energy Syst._ , vol. 45, no. 1, pp. 501–508, 2013.
  * [8] R. A. Jabr, “Polyhedral formulations and loop elimination constraints for distribution network expansion planning,” _IEEE Trans. Power Syst._ , vol. 28, no. 2, pp. 1888–1897, 2013.
  * [9] B. Kocuk, H. Jeon, S. S. Dey, J. Linderoth, J. Luedtke, and A. Sun, “A cycle-based formulation and valid inequalities for dc power transmission problems with switching,” _Oper. Res._ , vol. 64, no. 4, pp. 922 – 938, 2016\.
  * [10] A. Kyriacou, P. Demetriou, C. Panayiotou, and E. Kyriakides, “Controlled islanding solution for large-scale power systems,” _IEEE Trans. Power Syst._ , vol. 33, no. 2, pp. 1591–1602, 2017.
  * [11] G. Patsakis, D. Rajan, I. Aravena, and S. Oren, “Strong mixed-integer formulations for power system islanding and restoration,” _IEEE Trans. Power Syst._ , vol. 34, no. 6, pp. 4880–4888, 2019.
  * [12] L. Ding, F. M. Gonzalez-Longatt, P. Wall, and V. Terzija, “Two-step spectral clustering controlled islanding algorithm,” _IEEE Trans. Power Syst._ , vol. 28, no. 1, pp. 75–84, 2013.
  * [13] S. Fattahi, J. Lavaei, and A. Atamtürk, “A bound strengthening method for optimal transmission switching in power systems,” _IEEE Trans. Power Syst._ , vol. 34, no. 1, pp. 280–291, 2019.
  * [14] I. Tyuryukanov, M. A. M. M. van der Meijden, V. Terzija, and M. Popov, “Spectral mst-based graph outlier detection with application to clustering of power networks,” in _Proc. 2018 Power Systems Computation Conference (PSCC)_ , 2018, pp. 1 – 8.
  * [15] S. Chopra and C.-Y. Tsai, “Polyhedral approaches for the steiner tree problem on graphs,” in _Steiner Trees in Industry_ , ser. COOP, X. Z. Cheng and D.-Z. Du, Eds. Springer, 2001, vol. 11, pp. 175–201.
  * [16] A. Borghetti, “A mixed-integer linear programming approach for the computation of the minimum-losses radial configuration of electrical distribution networks,” _IEEE Trans. Power Syst._ , vol. 27, no. 3, pp. 1264–1273, 2012\.
  * [17] C. Hojny, I. Joormann, H. Lüthen, and M. Schmidt, “Mixed-integer programming techniques for the connected max-k-cut problem,” _Math. Program. Comput._ , vol. 13, pp. 75–132, 2021.
  * [18] Gurobi Optimization, LLC, “Gurobi optimizer reference manual,” 2020. [Online]. Available: http://www.gurobi.com
  * [19] S. Fliscounakis, P. Panciatici, F. Capitanescu, and L. Wehenkel, “Contingency ranking with respect to overloads in very large power systems taking into account uncertainty, preventive, and corrective actions,” _IEEE Trans. Power Syst._ , vol. 28, no. 4, pp. 4909–4917, 2013.
  * [20] I. Tyuryukanov, M. Popov, M. van der Meijden, and V. Terzija, “Slow coherency identification and power system dynamic model reduction by using orthogonal structure of electromechanical eigenvectors,” _IEEE Trans. Power Syst._ , vol. 36, no. 2, pp. 1482 – 1492, 2021.
  * [21] P. M. Anderson and A. A. Fouad, _Power system control and stability_ , 2nd ed., ser. IEEE Press power engineering series. Piscataway, N.J.: IEEE Press and Wiley-Interscience, 2003.

TABLE VI: Islands’ Power Imbalance Minimization using (6) $n$ | $K$ |  $UB$, | $g,~{}\%$ |  $P_{LS}$, |  $P_{\Delta}$, |  $P_{GS}$, |  $p_{ij}^{\Sigma}$,
---|---|---|---|---|---|---|---
|  | p.u. |  $T$, s | p.u. | p.u. | p.u. | p.u.
89 | 2 | 7.6364 | 19 s | 2.43 | 6.24 | 3.81 | 55.5
89 | 3 | 7.8258 | 52 s | 2.43 | 6.24 | 3.81 | 74.4
89 | 4 | 13.485 | 12 s | 4.88 | 11.1 | 6.26 | 66.6
89 | 5 | 13.488 | 12 s | 4.88 | 11.1 | 6.26 | 66.9
1354 | 2 | 17.023 | 0.4 s | -0.00 | 16.6 | 16.6 | 22.2
1354 | 3 | 17.38 | 11 s | -0.00 | 16.6 | 16.6 | 57.9
1354 | 4 | 17.123 | 1.0 s | -0.00 | 16.6 | 16.6 | 32.3
1354 | 5 | 17.604 | 76 s | -0.00 | 16.6 | 16.6 | 80.3
1354 | 6 | 18.089 | 180 s | -0.00 | 16.6 | 16.6 | 129
1354 | 7 | 21.668 | 17 % | 1.42 | 19.5 | 18.1 | 155
1354 | 8 | 18.565 | 350 s | -0.00 | 16.6 | 16.6 | 176
1888 | 2 | 9.9447 | 0.5 s | -0.00 | 9.81 | 9.81 | 3.97
1888 | 3 | 10.821 | 65 s | -0.00 | 9.81 | 9.81 | 91.7
1888 | 4 | 10.431 | 6.7 s | 0.03 | 9.87 | 9.84 | 45.1
1888 | 5 | 10.68 | 79 s | -0.00 | 9.81 | 9.81 | 77.5
1888 | 6 | 10.92 | 170 s | -0.00 | 9.81 | 9.81 | 101
1888 | 7 | 11.298 | 150 s | -0.00 | 9.81 | 9.81 | 139
1888 | 8 | 11.208 | 283 s | -0.00 | 9.81 | 9.81 | 131