# Detection and prediction of clopidogrel treatment failures
using longitudinal structured electronic health records

Samuel Kim, In Gu Sean Lee, Mijeong Irene Ban, and Jane Chiang

Cipherome Inc., San Jose, CA, U.S.A.
{sam.kim, sean.lee, irene.ban<EMAIL_ADDRESS>
###### Abstract

We propose machine learning algorithms to automatically detect and predict
clopidogrel treatment failure using longitudinal structured electronic health
records (EHR). By drawing analogies between natural language and structured
EHR, we introduce various machine learning algorithms used in natural language
processing (NLP) applications to build models for treatment failure detection
and prediction. In this regard, we generated a cohort of patients with
clopidogrel prescriptions from UK Biobank and annotated if the patients had
treatment failure events within one year of the first clopidogrel
prescription; out of 502,527 patients, 1,824 patients were identified as
treatment failure cases, and 6,859 patients were considered as control cases.
From the dataset, we gathered diagnoses, prescriptions, and procedure records
together per patient and organized them into visits with the same date to
build models. The models were built for two different tasks, i.e., detection
and prediction, and the experimental results showed that time series models
outperform bag-of-words approaches in both tasks. In particular, a
Transformer-based model, namely BERT, could reach 0.928 AUC in detection tasks
and 0.729 AUC in prediction tasks. BERT also showed competence over other time
series models when there is not enough training data, because it leverages the
pre-training procedure using large unlabeled data.

##  I Introduction

Adverse drug reactions (ADRs) and treatment failures (TF), both of which can
occur due to inadequate concentrations of medications in the system,
significantly impact patients’ lives and are serious burdens to the healthcare
system. ADRs, which may occur from excess active substance in the system, are
the fourth leading cause of death behind cardiovascular diseases. Treatment
failures may occur when therapeutic doses are not reached. Subsequent
hospitalizations may significantly increase healthcare costs [1]. The variable
presentation makes identifying drug responses difficult, and the impact is
likely underrepresented [2]. The multi-factorial nature of drug responses adds
to the complexity; patient’s age, renal or liver function, comorbidities,
comedications, lifestyle, and genetic predispositions are all important
factors affecting drug responses [3, 4].

Among multitude of drugs, we studied clopidogrel, a P2Y12 inhibitor, one of
the most widely prescribed drug in the US. Its use with aspirin is commonly
referred to as dual antiplatelet therapy (DAPT), and is the standard of care
following stent placement for cardiovascular disease. However, various drug
responses to this drug has been reported in several studies raising concerns
for its safety [5]. High concentration of active metabolite can cause
bleeding, while low concentration of active metabolite can cause recurrent
thrombosis, both of which can lead to medical emergencies. There is a clear
need for a screening tool to better predict which patients may be at risk for
such events.

Therefore, we focused on building machine learning models to automatically
detect and predict the treatment failures on clopidogrel. More specifically,
we used structured electronic health records (EHR) which contain information
on many of the risk factors aforementioned. Easier access to medical records
have enabled numerous applications using machine learning approaches. We used
only structured EHR, in particular data generated by clinical coding schemes.
Clinical episodes are coded by clinical experts to describe patients’ health
condition and treatments. Codes include disease diagnoses, procedures, and
medicine prescriptions along with the dates.

As structured EHR is inherently sequential, many researchers have applied time
series modeling strategies. For example, Choi et al. proposed a recurrent
neural network (RNN) architecture to detect heart failures [6]. The authors
applied an attention mechanism on grouped diagnostic codes in a reverse time
order and showed clinically interpretable results. Li et al. proposed BEHRT to
learn patterns of diagnostic codes and used the model to predict unseen
diagnoses [7]. Although the original codes were generated in Read codes [8]
and ICD-10 codes [9], they were mapped into Caliber codes [10] so that the
vocabulary size became 301. Shang et al. proposed a combined architecture of
Graph Neural Networks (GNNs) [11] and Bidirectional Encoder Representations
from Transformer (BERT) [12] for medication recommendation [13]. The authors
used the GNNs to represent ontology embeddings and the BERT to generate visit
embeddings for downstream applications. Ramsy et al. proposed Med-BERT to
build a disease classification such as heart failure among patients with
diabetes (DHF) and pancreatic cancer (PaCa) [14]. The model was trained using
ICD-9/10 codes as tokens directly and fined-tuned for disease classification.

In this paper, we investigated various time series modeling approaches to
model the longitudinal structured EHR toward treatment failures detection and
prediction. We hypothesized that the dynamic nature of medical records needs
to be modeled. In that regard, we annotated clopidogrel treatment failures
using publicly available data and built prediction models and detection
models. Although previous related works are somewhat limited by only using
diagnosis information, we proposed to use all possible data in addition to
diagnoses, i.e., procedures and medication prescriptions. We hypothesized that
they are complementary to represent patients’ medical condition and treatments
and modeling the interactions of them would provide better discriminating
power.

##  II Data

###  II-A Dataset

We used the UK Biobank which consists of data collected from 502,527
participants [15]. Volunteers aged 40 to 70 were recruited from England,
Scotland, and Wales and invited to assessment centers between 2006 and 2010.
Data collected during the visits include biosamples, physical examination
measurements, and questionnaire answers followed by interviews. Genomic data,
both sequencing and genotyping, have been generated using the biosamples
collected as well. Also, an extensive and comprehensive medical history data
have also been made available through hospital in patient and primary care
data from external data sources.

We extracted prescriptions, diagnoses and procedure records, along with dates,
for all participants from the UK Biobank. All prescription records were from
general practitioner (GP) data coded in Read and British National Formulary
(BNF) depending on the data supplier. Diagnoses and procedure records were
from hospital in-patient records only and coded in ICD-9/10 and OPCS-3/4
respectively.

###  II-B Annotation

Treatment failure (TF) was defined as having a TF event within one year of the
very first clopidogrel prescription. TF events include ischemic stroke,
myocardial infarct, stent thrombosis and recurrent thrombosis or stenting, all
of which were identified using ICD-9/10 and OPCS-3/4 codes. Clopidogrel
prescriptions were identified using the substance or brand names of
clopidogrel or respective read codes.

Subjects with events occurring within 7 days of the first prescription were
excluded as it is unclear whether those events were associated with
clopidogrel. The visit had to be through the emergency room to be valid in
order to exclude follow up visits from previous events. From the dataset, we
found 9,867 subjects with clopidogrel prescriptions. Among them, we labeled
1,824 patients as treatment failure cases and 6,859 as control cases; 1,184
subjects were excluded due to data inconsistencies or ambiguities.

##  III Methodologies

###  III-A Data Processing

We gathered various codes together per subject and organized them into visits
with the same date. As a result, a subject has multiple visits and a visit has
multiple codes. Table I shows the statistics of the codes and the visits.

TABLE I: Statistics of structured EHR data in UKBiobank. # of subjects | 502,527
---|---
# of subjects with 1+ visits | 465,506
# of unique visits | 22,063,417
# of unique codes | 31,589
Average # of visits per subject (max/med/min) | 58.3 (2805 / 13 / 1)
Average # of codes per visit (max/med/min) | 2.7 (61 / 2 / 1)

Note that the number of codes per visit is relatively small compared to the
number of words per sentence in natural language processing (NLP)
applications. Also, the codes within in a visit do not have obvious
ordinality. On the other hand, the number of visits per patient is fairly
large and their order is less ambiguous due to timestamps. Therefore, we
flattened the codes of the visits and concatenated them across all visits so
that a patient has a ordered sequence of codes as illustrated in Fig. 1.

Figure 1: Example of a patient’s longitudinal structured EHR on timeline and
how we process them.

###  III-B Model Architectures

We used various machine learning algorithms including gated recurrent neural
network (GRU) [16], and long short-term memory (LSTM) [17], and BERT [12] for
time series modeling. We also used random forest (RF) and logistic regression
(LR) with bag-of-words approaches as baselines.

For recurrent-based neural network approaches such GRU and LSTM, we
empirically chose one hidden layer of GRU or LSTM with 768 units. For BERT
models, we used the original architecture in [18] for simplicity. In both
approaches, only the most recent 512 codes per subject were used if there were
more. To build a pre-trained model for BERT, we used the data from the
subjects that were not annotated with treatment failures. Then, we fine-tuned
the pre-trained model using the labeled data. We replaced the top layer of the
pre-trained model with a fully connected layer with one node with binary
cross-entropy as a loss function.

We also built a tokenizer that converts individual clinical codes into unique
integer values to feed into machine learning algorithms. The tokenizer is
built using all the codes together unless noted otherwise.

###  III-C Experimental Setup

We designed two experimental tasks: 1) detecting treatment failures within a
year after the first prescription of clopidogrel and 2) predicting treatment
failures within a year after the first prescription of clopidogrel. The main
difference of these two tasks is the direction of point of views on the
timeline. The former is looking back from one year after the first
prescription if there have been any treatment failures, while the latter is
looking forward at the first prescription date to predict if there will be any
treatment failures. From the data perspectives, the former has an access to
the all data until one year after the first prescription date while the latter
only has access the data before the first prescription date (see Fig. 1 for an
example).

We randomly split the labeled data into training and testing data with 80/20
ratio. This split was used across all experiments for fair comparisons. Twenty
percent of the training data was used as validation data if necessary. We used
receiver operating characteristic (ROC) curves and their areas under curve
(AUC) to measure the performance of each experiment.

##  IV Results and Discussion

Fig. 2 (a) and (b) illustrates the performance of the models for the detection
and prediction tasks. As shown in the figures, the results of detection tasks
are better than the prediction tasks. This is reasonable because the detection
tasks are to make decisions with more data after the prescription of the drug
while the prediction tasks are to predict only with the data before the
prescription.

(a) Detection tasks

(b) Prediction tasks

Figure 2: Performance of the models for adverse drug reactions in terms of ROC
and its AUC: (a) detection tasks and (b) prediction tasks.

(a) Detection tasks

(b) Prediction tasks

Figure 3: Performance of the models for adverse drug reactions with respect to
number of training samples in terms of AUC.

It is notable that using time series modeling strategies such as BERT, LSTM
and GRU outperforms the ones using bag-of-words approaches such as RF and LR
in both tasks. This vindicates that the dynamics in the longitudinal medical
records convey significant information regarding drug responses. Note that the
best model for detection tasks is GRU (AUC=0.938) rather than BERT
(AUC=0.928). Although we initially expected BERT to be the best model as it
used a large unlabeled dataset to build the pre-trained model and fine-tuned
with the labeled data, it was somewhat difficult to observe benefits of the
pre-trained model.

TABLE II: Performance of the models for adverse drug reactions with different
types of data in terms of AUC. Checkmarks are the modalities that were used to
build the corresponding model.

(a) Detection tasks
---
Proc. | Diag. | Pres. | BERT | LSTM | GRU | RF | LR
|  |  | 0.928 | 0.838 | 0.938 | 0.798 | 0.792
|  |  | 0.716 | 0.805 | 0.801 | 0.779 | 0.793
|  |  | 0.759 | 0.810 | 0.809 | 0.771 | 0.769
|  |  | 0.721 | 0.723 | 0.722 | 0.725 | 0.628
|  |  | 0.829 | 0.820 | 0.834 | 0.810 | 0.813
|  |  | 0.877 | 0.822 | 0.882 | 0.794 | 0.776
|  |  | 0.808 | 0.726 | 0.887 | 0.771 | 0.750

(b) Prediction tasks
---
Proc. | Diag. | Pres. | BERT | LSTM | GRU | RF | LR
|  |  | 0.729 | 0.703 | 0.692 | 0.669 | 0.612
|  |  | 0.606 | 0.632 | 0.624 | 0.594 | 0.602
|  |  | 0.604 | 0.623 | 0.613 | 0.643 | 0.628
|  |  | 0.615 | 0.648 | 0.691 | 0.669 | 0.563
|  |  | 0.622 | 0.663 | 0.666 | 0.659 | 0.640
|  |  | 0.605 | 0.667 | 0.662 | 0.672 | 0.581
|  |  | 0.615 | 0.679 | 0.664 | 0.675 | 0.600

Therefore, we performed extra experiments with different training data sizes.
We randomly chose a subset of training data with variable sizes to simulate
the cases where there were not sufficient training data. Note that we kept the
same test dataset across the experiments for fair comparison. Fig. 3 depicts
AUCs of the models with respect to number of training samples. As shown in the
figures, it is evident that the AUCs of sequential modeling approaches (solid
lies) degrade as the number of training samples decreases while bag-of-words
approaches (dotted lines) do not have clear patterns. This is partially
because the number of parameters of sequential models are far greater than the
ones of bag-of-words models so they suffer from lack of data. It is also
notable that BERT outperforms other models with smaller training datasets
particularly in detection tasks. We argue that BERT benefits from leveraging
the pre-training procedure and quickly adopting their parameters with the
labeled training data.

We also performed ablation experiments which enable us to analyze the impact
of different types of medical records. For these experiments, we built models
for individual modalities, i.e. procedures, diagnoses, and prescriptions, and
also combinations. We also built separate tokenizers for each experiment. As
shown in Table II, the models with all the modalities outperform the ones with
partial information especially using sequential models. This confirms that
modeling the interactions in multiple modalities in structured EHR provides
better predictive power.

##  V Conclusions

We introduced ways to model longitudinal structured EHR toward detecting and
predicting treatment failures on clopidogrel. We applied time series modeling
strategies that are commonly used in NLP applications. The experimental
results showed that the algorithms that model time dynamics outperformed bag-
of-words approaches and suggested to use pre-trained models when the amount of
labeled data is small. It was also shown that using various types of medical
records together improved the performance in detecting and predicting
treatment failures.

In future, we will investigate visit embeddings to represent all the codes
within a visit to mitigate ordinality issues. We also want to incorporate with
lab data such as glomerular filtration rate, liver function test, and
hemoglobin levels, to improve model performance.

##  VI Acknowledgement

This research has been conducted using the UK Biobank Resource under
Application Number 52031.

## References

  * [1] Ahmedin Jemal, Elizabeth Ward, Yongping Hao, and Michael Thun,  “Trends in the Leading Causes of Death in the United States, 1970-2002,”  JAMA, vol. 294, no. 10, 09 2005.
  * [2] Lorna Hazell and Saad A.W. Shakir,  “Under-reporting of adverse drug reactions,”  Drug Safety, vol. 29, no. 5, pp. 385–396, 2006.
  * [3] Steven M. Handler, Joseph T. Hanlon, Subashan Perera, et al.,  “Consensus list of signals to detect potential adverse drug reactions in nursing homes.,”  Journal of the American Geriatrics Society, vol. 56, no. 5, pp. 808–815, 2008.
  * [4] U. A. Meyer,  “Pharmacogenetics and adverse drug reactions.,”  Lancet (London, England), vol. 356, no. 9242, pp. 1667–1671, 2000\.
  * [5] M. Cattaneo,  “Resistance to antiplatelet drugs: molecular mechanisms and laboratory detection,”  Journal of Thrombosis and Haemostasis, vol. 5, no. s1, pp. 230–237, 2007.
  * [6] Edward Choi, Mohammad Taha Bahadori, Joshua A. Kulas, Andy Schuetz, Walter F. Stewart, and Jimeng Sun,  “RETAIN: An Interpretable Predictive Model for Healthcare using Reverse Time Attention Mechanism,”  arXiv:1608.05745 [cs], Feb. 2017.
  * [7] Laila Rasmy, Yang Xiang, Ziqian Xie, Cui Tao, and Degui Zhi,  “Med-BERT: pretrained contextualized embeddings on large-scale structured electronic health records for disease prediction,”  npj Digital Medicine, vol. 4, no. 1, pp. 86, Dec. 2021.
  * [8] National Health Service,  “Read codes,” https://digital.nhs.uk/services/terminology-and-classifications/read-codes.
  * [9] World Health Organization,  “ICD-10 : international statistical classification of diseases and related health problems: tenth revision,” 2004.
  * [10] Valerie Kuan et al.,  “A chronological map of 308 physical and mental health conditions from 4 million individuals in the english national health service,”  The Lancet Digital Health, vol. 1, 05 2019.
  * [11] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and Philip S. Yu,  “A Comprehensive Survey on Graph Neural Networks,”  IEEE Transactions on Neural Networks and Learning Systems, vol. 32, no. 1, pp. 4–24, Jan. 2021,  arXiv: 1901.00596.
  * [12] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova,  “BERT: pre-training of deep bidirectional transformers for language understanding,”  arXiv:1810.04805 [cs], May 2019.
  * [13] Junyuan Shang, Tengfei Ma, Cao Xiao, and Jimeng Sun,  “Pre-training of Graph Augmented Transformers for Medication Recommendation,”  arXiv:1906.00346 [cs], Nov. 2019.
  * [14] Yikuan Li, Shishir Rao, José Roberto Ayala Solares, et al.,  “BEHRT: Transformer for electronic health records,”  Scientific Reports, vol. 10, no. 1, pp. 7155, 2020.
  * [15] Cathie Sudlow et al.,  “Uk biobank: An open access resource for identifying the causes of a wide range of complex diseases of middle and old age,”  PLOS Medicine, vol. 12, no. 3, pp. 1–10, 03 2015.
  * [16] Junyoung Chung, Caglar Gulcehre, Kyunghyun Cho, and Yoshua Bengio,  “Empirical evaluation of gated recurrent neural networks on sequence modeling,”  in NIPS 2014 Workshop on Deep Learning, December 2014, 2014.
  * [17] Sepp Hochreiter and Jürgen Schmidhuber,  “Long short-term memory,”  Neural computation, vol. 9, pp. 1735–80, 12 1997.
  * [18] Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, et al.,  “Transformers: State-of-the-art natural language processing,”  in Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations, Online, Oct. 2020, pp. 38–45, Association for Computational Linguistics.


# Positive solutions for non-variational fractional elliptic systems with
negative exponents 1112000 Mathematics Subject Classification: 35R11; 35B25;
35A16. 222Key words: Fractional Laplace operator, non-variational elliptic
systems, negative exponents.

Anderson L. A. de Araujo 333E-mail addresses: anderson.araujo@ufv.br(A.
Araujo)444A.L.A de Araujo was partially supported by FAPEMIG/FORTIS.
Departamento de Matemática, Universidade Federal de Viçosa,
CCE, 36570-000, Viçosa, MG, Brazil
Luiz F. O. Faria 555E-mail addresses: luiz.faria@ufjf.edu.br(L. F. O.
Faria)666L. F.O Faria was partially supported by FAPEMIG CEX APQ 02374/17.
Departamento de Matemática, Universidade Federal de Juiz de Fora,
ICE, 36036-330, Juiz de Fora, MG, Brazil
Edir Junior F. Leite 777E-mail addresses<EMAIL_ADDRESS>(E.J.F. Leite)
Departamento de Matemática, Universidade Federal de Viçosa,
CCE, 36570-000, Viçosa, MG, Brazil
Olímpio H. Miyagaki 888E-mail addresses<EMAIL_ADDRESS>(O. H.
Miyagaki) 999O. H. M. was supported in part by INCTMAT/BRAZIL and CNPQ/BRAZIL
PROC. 307061/2018-3
Departamento de Matemática, Universidade Federal de Juiz de Fora,
ICE, 36036-330, Juiz de Fora, MG, Brazil

Abstract

In this paper, we study strongly coupled elliptic systems in non-variational
form with negative exponents involving fractional Laplace operators. We
investigate the existence, nonexistence, and uniqueness of the positive
classical solution. The results obtained here are a natural extension of the
results obtained by Ghergu, in [6], for the fractional case.

##  1 Introduction and main results

The present paper deals with existence, nonexistence, and uniqueness of
positive solutions for elliptic systems of the form

$\left\\{\begin{array}[]{llll}(-\Delta)^{s}u=u^{-p}v^{-q}&{\rm in}\ \
\Omega\\\ (-\Delta)^{t}v=u^{-r}v^{-\theta}&{\rm in}\ \ \Omega\\\ u=v=0&{\rm
in}\ \ \mathbb{R}^{n}\setminus\Omega\end{array}\right.$ (1)

where $\Omega$ is a smooth bounded open subset of $\mathbb{R}^{n}$, $n\geq 2$,
$0<s,t<1$, $r,q>0$, $p,\theta\geq 0$ and the fractional Laplace operator
$(-\Delta)^{s}$ is defined as

$(-\Delta)^{s}u(x)=C(n,s)\lim_{\varepsilon\searrow
0}\int\limits_{\mathbb{R}^{n}\setminus
B_{\varepsilon}(x)}\frac{u(x)-u(y)}{|x-y|^{n+2s}}\;dy\,,$

for all $x\in\mathbb{R}^{n}$ and

$C(n,s)=\left(\int\limits_{\mathbb{R}^{n}}\frac{1-\cos(\zeta_{1})}{|\zeta|^{n+2s}}\;d\zeta\right)^{-1}$

with $\zeta=(\zeta_{1},\ldots,\zeta_{n})\in\mathbb{R}^{n}$. A natural space
for this operator is a weighted $L_{1}$-space:

$L_{s}:=\left\\{u:\mathbb{R}^{n}\rightarrow\mathbb{R}:\int\limits_{\mathbb{R}^{n}}\frac{|u(x)|}{1+|x|^{n+2s}}dx<+\infty\right\\}.$

The norm in $L_{s}$ is naturally given by

$\|u\|_{L_{s}}=\int\limits_{\mathbb{R}^{n}}\frac{|u(x)|}{1+|x|^{n+2s}}dx.$

The study of system (1) was mainly motivated from the well known fractional
Lane-Emden problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&u^{p}&{\rm in}\ \ \Omega\\\
u&=&0&{\rm in}\ \ \mathbb{R}^{n}\setminus\Omega\end{array}\right.,$ (2)

where $\Omega$ is a smooth bounded open subset of $\mathbb{R}^{n}$, $n\geq 1$
and $0<s<1$.

Recently, it has been proved in [16] that this problem admits at least one
positive solution for $1<p<\frac{n+2s}{n-2s}$. The nonexistence has been
established in [14] whenever $p\geq\frac{n+2s}{n-2s}$ and $\Omega$ is star-
shaped. These results were known long before for $s=1$, see the classical
references [4, 5, 10].

For system of the type (1) with $p=0=\theta$ and $r,q<0$, existence results of
positive solutions have been established when $qr>1$ in [7] for $s\neq t$ and
in [8] for $s=t$. The latter also proves existence and uniqueness of positive
solution in the case that $qr<1$. Finally, when $qr=1$, the behavior of (1) is
resonant and the related eigenvalue problem has been studied in [9].

Nowadays, there has been some interest in systems of the type (1) with
$p,\theta\geq 0$ and $q,r>0$. In [6], the author studied existence,
nonexistence, uniqueness, and regularity of solutions for the system (1) with
$s=1=t$.

In this paper, we are going to treat the system (1) in the case $p,\theta\geq
0$ and $q,r>0$. In this structure, the system above corresponds to the
prototype equation (2) in which the exponent $p$ is negative and generalize
the results obtained in [6]. It is well known that for such a range of
exponents, the system (1) does not have a variational structure. To overcome
this, we employ the sub-super method, which our approach relies on the
boundary behavior of solutions to (2) (with $p<0$) or more generally, to
singular elliptic problems of the type

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&K(x)u^{-p},\ \ u>0&{\rm in}\ \
\Omega\\\ u&=&0&{\rm in}\ \ \mathbb{R}^{n}\setminus\Omega\end{array}\right.$
(3)

where $K\in C^{\nu}_{loc}(\Omega),\nu\in(0,1)$, such that
$\inf\limits_{\Omega}K>0$ and satisfies for some $0\leq q<2s$ and
$C_{1},C_{2}>0$

$C_{1}d(x)^{-q}\leq K(x)\leq C_{2}d(x)^{-q}\text{ in }\Omega,$

where $d(x)=dist(x,\partial\Omega)$, studied by Adimurthi, Giacomoni and
Santra in [3].

We say that a pair $(u,v)$ of continuous function in $\Omega$ and bounded in
$\mathbb{R}^{n}$ is a positive classical solution of system (1), if
$(-\Delta)^{s}u(x)$ and $(-\Delta)^{t}v(x)$ are well defined for all
$x\in\Omega$, further $u$ and $v$ are positive in $\Omega$ and all equalities
in (1) hold pointwise in each corresponding set. Positive classical super and
subsolutions are defined similarly.

We will establish our first result concerning the system (1).

######  Theorem 1.1.

(Nonexistence). Let $p,\theta\geq 0$, $r,q>0.$ Then the system (1) has no
positive classical solutions, provided that one of the following conditions
holds:

  * $(i)$

$\frac{qt}{s}+p<1$ and $r\geq\frac{2t}{s}$;

  * $(ii)$

$\frac{qt}{s}+p>1$ and $r(2s-qt)\geq 2s(1+p)$;

  * $(iii)$

$\frac{rs}{t}+\theta<1$ and $q\geq\frac{2s}{t}$;

  * $(iv)$

$\frac{rs}{t}+\theta>1$ and $q(2t-rs)\geq 2t(1+\theta)$;

  * $(v)$

$p>\max\\{1,\frac{rs}{t}-1\\},\frac{2rs}{t}>(1-\theta)(1+p)$ and
$qt(1+p-\frac{rs}{t})>(1+p)(1+\theta)s$;

  * $(vi)$

$\theta>\max\\{1,\frac{qt}{s}-1\\},\frac{2qt}{s}>(1-p)(1+\theta)$ and
$rs(1+\theta-\frac{qt}{s})>(1+p)(1+\theta)t$.

######  Remark 1.1.

The conditions (i) and (ii) in Theorem 1.1 impose conditions on the exponent
$q$ to vary on the interval $(0,\frac{2s}{t})$, while in (v) the exponent $q$
can take any value greater than $\frac{2s}{t}$, provided adjusting the other
three exponents $p,r,\theta$ conveniently. Finally, from the conditions (iii),
(iv) and (vi), the exponent $r$ is also restricted as above.

Define the following quantities

$\alpha=p+\frac{qt}{s}\min\left\\{1,\frac{2t-rs}{(1+\theta)t}\right\\},\ \ \
\beta=\theta+\frac{rs}{t}\min\left\\{1,\frac{2s-qt}{(1+p)s}\right\\}.$

These above quantities $\alpha$ and $\beta$ are related to the boundary
behavior of the solution to the singular elliptic problem (3), as they will be
explained in Proposition 2.3 below.

Next, we will state the existence of classical solutions to (1).

######  Theorem 1.2.

(Existence). Let $p,\theta\geq 0$, $q,r>0$ satisfying the inequality

$(1+p)(1+\theta)-qr>0.$ (4)

In addition, assume that one of the following conditions below holds:

  * $(i)$

$\alpha\leq 1$ and $r<\frac{2t}{s}$;

  * $(ii)$

$\beta\leq 1$ and $q<\frac{2s}{t}$;

  * $(iii)$

$p,\theta\geq 1$, $r<\frac{2t}{s}$ and $q<\frac{2s}{t}$.

Then, the system (1) has at least one positive classical solution
$(u,v)\in(C^{\eta}(\mathbb{R}^{n}))^{2}$, for some $\eta\in(0,1)$.

The proof is made invoking the Schauder’s fixed point theorem in a suitable
chosen closed convex subset of $(C^{\eta}(\mathbb{R}^{n}))^{2}$, for some
$\eta\in(0,1)$, which contains all the functions having a certain rate of
decay expressed in terms of the distance function $d(x)$ up to the boundary of
$\Omega$.

The following necessary and sufficient conditions for the existence of
classical solutions to (1) follows directly from Theorem 1.1(i) and (iii) and
Theorem 1.2(i) and (ii).

######  Corollary 1.1.

Let $p,\theta\geq 0$, $q,r>0$ satisfy (4).

  * $(i)$

Assume $\frac{qt}{s}+p<1$. Then system (1) has positive classical solutions if
and only if $r<\frac{2t}{s}$;

  * $(ii)$

Assume $\frac{rs}{t}+\theta<1$. Then system (1) has positive classical
solutions if and only if $q<\frac{2s}{t}$.

######  Theorem 1.3.

(Uniqueness). Let $p,\theta\geq 0$, $q,r>0$, satisfy (4) and one of the
following conditions:

  * $(i)$

$\frac{qt}{s}+p<1$ and $r<\frac{2t}{s}$;

  * $(ii)$

$\frac{rs}{t}+\theta<1$ and $q<\frac{2s}{t}$.

Then, the system (1) has a unique positive classical solution.

Several methods have been employed in the proof of existence, nonexistence and
uniqueness results of positive solutions of elliptic systems. Our approach is
inspired by a method developed by Ghergu in [6] to treat systems involving
Laplace operators based on boundary behavior of the solution to (3), when
$s=1$. Particularly, the boundary behavior of the solution to (3), proved by
Adimurthi, Giacomoni, and Santra [3], as well as some fundamental results to
be proved in the next section will play an important role in the proofs of
Theorems of this work.

The paper is organized as follows. In Section 2 we obtain some preliminary
properties related to the boundary behavior of the solution to (3). The rest
of the Sections are devoted to the proofs of our results.

##  2 Notation and auxiliary results

Consider the nonlocal eigenvalue problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&\lambda u&{\rm in}\ \ \Omega\\\
u&=&0&{\rm in}\ \ \mathbb{R}^{n}\setminus\Omega\end{array}\right..$ (5)

Since the operator $(-\Delta)^{s}$ is self-adjoint, by using a weak
formulation and a suitable variational framework, Servadei and Valdinoci [15]
investigated in detail the discrete spectrum of $(-\Delta)^{s}$ in $\Omega$
for any $s\in(0,1)$. In particular, they proved that the first eigenvalue
$\lambda_{1}(s):=\lambda_{1}((-\Delta)^{s})$ is positive, simple and
characterized by

$\lambda_{1}(s)=\inf_{u\in
X(\Omega)\setminus\\{0\\}}\frac{\int\limits_{\mathbb{R}^{n}}|(-\Delta)^{\frac{s}{2}}u|^{2}dx}{\int\limits_{\mathbb{R}^{n}}|u|^{2}}dx,$

where

$X(\Omega):=\\{u\in H^{s}(\mathbb{R}^{n}):u=0\text{ a.e in
}\mathbb{R}^{n}\setminus\Omega\\}.$

Let $\varphi_{s}$ be a nonnegative eigenfunction corresponding to
$\lambda_{1}(s)$ in the weak sense. Results of Hölder regularity to the
operator $(-\Delta)^{s}$ obtained by Ros-Oton and Serra [12] imply that
$\varphi_{s}\in C^{s}(\mathbb{R}^{n})$ and moreover is a classical solution of
(5) which is positive in $\Omega$. The last claim follows from Silvestre’s
strong maximum principle [17] which holds for classical supersolutions
(subsolutions).

By suitable normalization we may assume $|\varphi_{s}|_{\infty}=1$. In
addition, it follows from the results in [13] that

$cd(x)^{s}\leq\varphi_{s}(x)\leq\frac{1}{c}d(x)^{s},$ (6)

for some positive constant $c$.

We denote by $G_{s}(\cdot,\cdot)$ the Green’s function of the fractional
Laplace operator $(-\Delta)^{s}$ on $\Omega$. Let $w$ be a weak solution of
the following problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}w&=&h&{\rm in}\ \ \Omega\\\
w&=&0&{\rm in}\ \ \mathbb{R}^{n}\setminus\Omega\end{array}\right..$ (7)

If $h\in C^{\alpha}_{loc}(\Omega)$, for some $\alpha\in(0,1)$, by Theorem 2.5
of [11], there exists $\gamma>0$ such that $w\in C^{2s+\gamma}_{loc}(\Omega)$
is a classical solution of (7), i.e, both equalities hold pointwise in each
corresponding set. Therefore,

$w(x)=\int\limits_{\Omega}G_{s}(x,y)h(y)dy\text{ in }\Omega\text{ and
}w(x)=0\text{ in }\mathbb{R}^{n}\setminus\Omega.$ (8)

Reciprocally, if $h\in C^{\alpha}(\overline{\Omega})$, for some
$\alpha\in(0,1)$, by Theorem 1.2.3 of [2] the function defined by setting (8)
belongs to $C^{2s+\varepsilon}_{loc}(\Omega)\cap C(\overline{\Omega})\cap
L_{s}$, fulfills $d(x)^{1-s}w\in C(\overline{\Omega})$, and $w$ is the only
classical solution of problem (7).

Now, let $\phi_{s}$ be the function that satisfies

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}\phi_{s}&=&1&{\rm in}\ \ \Omega\\\
\phi_{s}&=&0&{\rm in}\ \ \mathbb{R}^{n}\setminus\Omega\end{array}\right..$

By Silvestre’s strong maximum principle (see [17]), we get $\phi_{s}(x)>0$ in
$\Omega$. Therefore,

$\varphi_{s}(x)=\lambda_{1}(s)\int\limits_{\Omega}G_{s}(x,y)\varphi_{s}(y)dy$

and

$\phi_{s}(x)=\int\limits_{\Omega}G_{s}(x,y)dy,$

which, as a consequence of the normalization of $\varphi_{s}$, leads to

$\varphi_{s}\leq\lambda_{1}(s)\phi_{s}.$ (9)

An important tool for the uniqueness result of solutions of the system (1) is
as follows:

######  Proposition 2.1.

Let $p\geq 0$ and $\psi:\Omega\rightarrow(0,\infty)$ be a continuous function.
If $\underline{u}$ is a positive classical subsolution and $\overline{u}$ is a
positive classical supersolution of

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&\psi(x)u^{-p}&{\rm in}\ \
\Omega\\\ u&=&0&{\rm in}\ \ \mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

then $\underline{u}\leq\overline{u}$ in $\Omega$.

Proof. If $p=0$ the result is a consequence of the Silvestre’s strong maximum
principle. Suppose $p>0,$ and assume by contradiction that the set
$\omega:=\\{x\in\Omega:\overline{u}(x)<\underline{u}(x)\\}$ is not empty and
let $w:=\underline{u}-\overline{u}$. Then, $w$ achieves its maximum on
$\Omega$ at a point $x_{0}\in\omega$. Then,

$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle
C(n,s)\lim_{\varepsilon\searrow 0}\int\limits_{\mathbb{R}^{n}\setminus
B_{\varepsilon}(x_{0})}\frac{w(x_{0})-w(y)}{|x_{0}-y|^{n+2s}}\;dy=(-\Delta)^{s}w(x_{0})$
$\displaystyle\leq$
$\displaystyle\psi(x_{0})[\underline{u}(x_{0})^{-p}-\overline{u}(x_{0})^{-p}]<0,$

which is a contradiction. Therefore, $\omega=\emptyset$, that is,
$\underline{u}\leq\overline{u}$ in $\Omega$.

Now an important tool for the nonexistence and uniqueness results of solutions
of the system (1) is as follows:

######  Proposition 2.2.

Let $(u,v)$ be a positive classical solution of system (1). Then, there exists
a constant $c>0$ such that

$u(x)\geq cd(x)^{s}\text{ and }v(x)\geq cd(x)^{t}\text{ in }\Omega.$ (10)

Proof. Let $(u,v)$ be a positive classical solution of (1). By inequalities
(6) and (9), there is a constant $c_{0}>0$ such that $\phi_{s}(x)\geq
c_{0}d(x)^{s}$ and $\phi_{t}(x)\geq c_{0}d(x)^{t}$ in $\Omega$. Notice that
$(-\Delta)^{s}u\geq C=(-\Delta)^{s}(C\phi_{s})$ in $\Omega$, where
$C=\min\limits_{\Omega}\\{u^{-p}v^{-q}\\}>0$. Then, by Silvestre’s strong
maximum principle, we deduce $u(x)\geq C\phi_{s}(x)\geq cd(x)^{s}$ in $\Omega$
and similarly $v(x)\geq cd(x)^{t}$ in $\Omega$, where $c>0$ is a positive
constant.

The following result is a direct consequence of Silvestre’s strong maximum
principle, inequality (6) and Theorem 1.2 of [3]. This is the key tool for the
existence, nonexistence and uniqueness results of solutions of the system (1).

######  Proposition 2.3.

Let $p\geq 0$ and $\gamma>0$. There are constants $c,C>0$ such that any
positive classical subsolution $\underline{u}$ and any positive classical
supersolution $\overline{u}$ of problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&d(x)^{-\gamma}u^{-p}&{\rm in}\
\ \Omega\\\ u&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right.,$ (11)

satisfies:

  * $(i)$

$\underline{u}(x)\leq Cd(x)^{s}$ and $\overline{u}\geq cd(x)^{s}$ in $\Omega$,
if $\frac{\gamma}{s}+p<1$;

  * $(ii)$

$\underline{u}(x)\leq
Cd(x)^{s}\left(\ln\left(\frac{2}{\varphi_{s}}\right)\right)^{\frac{1}{1+p}}$
and $\overline{u}\geq
cd(x)^{s}\left(\ln\left(\frac{2}{\varphi_{s}}\right)\right)^{\frac{1}{1+p}}$
in $\Omega$, if $\frac{\gamma}{s}+p=1$;

  * $(iii)$

$\underline{u}(x)\leq Cd(x)^{\frac{2s-\gamma}{1+p}}$ and $\overline{u}\geq
cd(x)^{\frac{2s-\gamma}{1+p}}$ in $\Omega$, if $\frac{\gamma}{s}+p>1$ with
$0<\gamma<2s$.

Finally, Theorem 1.2(iii) of [3] also guarantees that the problem (11) has no
positive classical solution, if $\gamma\geq 2s$. Such claimed is important for
the proof of nonexistence results of positive classical solutions of the
system (1).

##  3 Proof of Theorem 1.1

Notice that the system (1) is invariant under the transform
$(u,v,p,q,r,\theta,s)\rightarrow(v,u,\theta,r,q,p,t)$, so that, we need to
prove only the cases (i), (ii) and (v).

Suppose that there exists $(u,v)$ a positive classical solution of system (1).
By Proposition 2.2, we can find $c>0$ such that (10) holds.

(i) $\frac{qt}{s}+p<1$ and $r\geq\frac{2t}{s}$. Using the estimate (10) in the
first equation of the system (1) we have

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&v^{-q}u^{-p}\leq
c_{1}d(x)^{-qt}u^{-p},\ \ \ u>0&{\rm in}\ \ \Omega\\\ u&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right.,$ (12)

for some $c_{1}>0$. By Proposition 2.3(i) we conclude $u(x)\leq c_{2}d(x)^{s}$
in $\Omega$, for some $c_{2}>0$. From this and (10), we have there exists
$c_{0},c_{3}>0$ such that $c_{0}d(x)^{-rs}\leq u^{-r}\leq c_{3}d(x)^{-rs}$ in
$\Omega$. Using the second equation of (1) we find

$\left\\{\begin{array}[]{rrll}(-\Delta)^{t}v&=&u^{-r}v^{-\theta},\ \ \
v>0&{\rm in}\ \ \Omega\\\ v&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right..$ (13)

According to Theorem 1.2(iii) of [3], this is impossible, since $rs\geq 2t$.

(ii) $\frac{qt}{s}+p>1$ and $r(2s-qt)\geq 2t(1+p)$. In the same manner as
above, $u$ satisfies the problem (12). From Proposition 2.3(iii), we now
deduce

$u(x)\leq c_{1}d(x)^{\frac{2s-qt}{1+p}}$

in $\Omega$, for some $c_{1}>0$. Since $\frac{qt}{s}+p>1$, we have
$\frac{2s-qt}{1+p}>s$. From this and (10) we deduce

$u(x)\geq cd(x)^{\frac{2s-qt}{1+p}}.$

Then, there are $c_{2},c_{3}>0$ such that

$c_{2}d(x)^{-\frac{r(2s-qt)}{1+p}}\leq u^{-r}\leq
c_{3}d(x)^{-\frac{r(2s-qt)}{1+p}}\text{ in }\Omega.$

Now, using the second equation of (1) we have $v$ is a classical solution of
problem (13), which is impossible in view of Theorem 1.2(iii) of [3], since
$\frac{r(2s-qt)}{1+p}\geq 2t$.

(v) Let $M=\sup\limits_{x\in\overline{\Omega}}v$. From the first equation of
the system (1) we find

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&\geq&M^{-q}u^{-p},\ \ \ u>0&{\rm
in}\ \ \Omega\\\ u&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right..$

From Proposition 2.3(iii), we have $u(x)\geq c_{1}d(x)^{\frac{2s}{1+p}}$ in
$\Omega$, for some $c_{1}>0$. Combining this estimate with the second equation
of (1) we have

$\left\\{\begin{array}[]{rrll}(-\Delta)^{t}v&\leq&c_{2}d(x)^{-\frac{2rs}{1+p}}v^{-\theta},\
\ \ v>0&{\rm in}\ \ \Omega\\\ v&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

where $c_{2}>0$. Since $\frac{2rs}{t(1+p)}+\theta>1$, again by Proposition
2.3(iii) we obtain that the function $v$ satisfies

$v(x)\leq c_{3}d(x)^{\frac{2t(1+p)-2rs}{(1+p)(1+\theta)}}\text{ in }\Omega,$

for some $c_{3}>0$. Since $\frac{2rs}{t}>(1-\theta)(1+p)$, we have
$\frac{2t(1+p)-2rs}{(1+p)(1+\theta)}>t$. From this and (10) we deduce

$u(x)\geq cd(x)^{\frac{2t(1+p)-2rs}{(1+p)(1+\theta)}}.$

Then, there exists $c_{4},c_{5}>0$ such that

$c_{4}d(x)^{-\frac{q(2t(1+p)-2rs)}{(1+p)(1+\theta)}}\leq v^{-q}\leq
c_{5}d(x)^{-\frac{q(2t(1+p)-2rs)}{(1+p)(1+\theta)}}.$

Now, using the first equation of (1) we have $u$ is a classical solution of
problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&v^{-q}u^{-p},\ \ \ u>0&{\rm
in}\ \ \Omega\\\ u&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

which contradicts Theorem 1.2(iii) of [3], since
$qt(1+p-\frac{rs}{t})>(1+p)(1+\theta)s$. Thus, the system (1) has no positive
classical solutions. This completes the proof of Theorem 1.1.

##  4 Proof of Theorem 1.2

(i) The proof is made in six cases according to bounded behavior of singular
elliptic problems of the type (3), as it was pointed out in Proposition 2.3.

Case 1: $\frac{rs}{t}+\theta>1$ and $\alpha=\frac{q(2t-rs)}{s(1+\theta)}+p<1$.
From Proposition 2.3(i) and (iii) there exist $0<c_{1}<1<c_{2}$ such that:

$\bullet$ Any positive classical subsolution $\underline{u}$ and any positive
classical supersolution $\overline{u}$ of the problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&d(x)^{-\frac{q(2t-rs)}{(1+\theta)}}u^{-p}&{\rm
in}\ \ \Omega\\\ u&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right.,$ (14)

satisfy

$\overline{u}(x)\geq c_{1}d(x)^{s}\text{ and }\underline{u}(x)\leq
c_{2}d(x)^{s}\text{ in }\Omega.\\\ $ (15)

$\bullet$ Any positive classical subsolution $\underline{v}$ and any positive
classical supersolution $\overline{v}$ of the problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{t}v&=&d(x)^{-rs}v^{-\theta}&{\rm in}\
\ \Omega\\\ v&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right.,$ (16)

satisfy

$\overline{v}(x)\geq c_{1}d(x)^{\frac{2t-rs}{1+\theta}}\text{ and
}\underline{v}(x)\leq c_{2}d(x)^{\frac{2t-rs}{1+\theta}}\text{ in }\Omega.$

We fix $0<m_{1}<1<M_{1}$ and $0<m_{2}<1<M_{2}$ such that

$M_{1}^{\frac{r}{1+\theta}}m_{2}\leq c_{1}<c_{2}\leq
M_{1}m_{2}^{\frac{q}{1+p}}$ (17)

and

$M_{2}^{\frac{q}{1+p}}m_{1}\leq c_{1}<c_{2}\leq
M_{2}m_{1}^{\frac{r}{1+\theta}}.$ (18)

Note that the above choice of $m_{i},M_{i}$ $(i=1,2)$ is possible in view of
(4).

Let $\varepsilon_{1}>0$ small enough. Here $X$ stands for the Banach space

$\\{(u,v)\in C^{s-\varepsilon_{1}}(\mathbb{R}^{n})\times
C^{\frac{2t-rs}{1+\theta}-\varepsilon_{1}}(\mathbb{R}^{n}):u=v=0\ {\rm in}\
\mathbb{R}^{n}\setminus\Omega\\}$

endowed with the product norm

$\|(u,v)\|_{X}:=\|u\|_{C^{s-\varepsilon_{1}}(\mathbb{R}^{n})}+\|v\|_{C^{\frac{2t-rs}{1+\theta}-\varepsilon_{1}}(\mathbb{R}^{n})}\,.$

Set

$\displaystyle\mathcal{A}:=\left\\{(u,v)\in
X:\begin{array}[]{c}m_{1}d(x)^{s}\leq u\leq M_{1}d(x)^{s}\\\ \text{ and }\\\
m_{2}d(x)^{\frac{2t-rs}{1+\theta}}\leq v\leq
M_{2}d(x)^{\frac{2t-rs}{1+\theta}}\text{ in }\Omega\end{array}\right\\}.$

For any $(u,v)\in\mathcal{A}$, let $(Tu,Tv)$ be the unique positive classical
solution of the decoupled system

$\left\\{\begin{array}[]{llll}(-\Delta)^{s}(Tu)=v^{-q}(Tu)^{-p}&{\rm in}\ \
\Omega\\\ (-\Delta)^{t}(Tv)=u^{-r}(Tv)^{-\theta}&{\rm in}\ \ \Omega\\\
Tu=Tv=0&{\rm in}\ \ \mathbb{R}^{n}\setminus\Omega\end{array}\right.$ (20)

and define

$\mathcal{F}:\mathcal{A}\rightarrow X\text{ by }\mathcal{F}(u,v)=(Tu,Tv)\text{
for any }(u,v)\in\mathcal{A}.$ (21)

It is proved in [3], the existence of positive classical solution $Tu\in
C^{s}(\mathbb{R}^{n})$ and $Tv\in C^{\frac{2t-rs}{1+\theta}}(\mathbb{R}^{n}),$
and the uniqueness of the positive weak solution in each equation of the
system (20). We define the space $X$ as subspace of

$C^{s-\varepsilon_{1}}(\mathbb{R}^{n})\times
C^{\frac{2t-rs}{1+\theta}-\varepsilon_{1}}(\mathbb{R}^{n}),$

for some $\varepsilon_{1}>0$ small enough, to ensure the compactness of the
operator $\mathcal{F}$ (see Step 2 below).

Therefore, if $\mathcal{F}$ has a fixed point in $\mathcal{A}$, then the
existence of a positive classical solution to system (1) follows. To this end,
we shall prove that $\mathcal{F}$ satisfies the conditions:

$\mathcal{F}(\mathcal{A})\subseteq\mathcal{A},\ \ \ \mathcal{F}\text{ is
compact and continuous}.$

Hence, by Schauder’s fixed point theorem we deduce that $\mathcal{F}$ has a
fixed point in $\mathcal{A}$, which is a positive classical solution to (1).

Step 1: $\mathcal{F}(\mathcal{A})\subseteq\mathcal{A}$. Take
$(u,v)\in\mathcal{A}$. From the inequality

$v\leq M_{2}d(x)^{\frac{2t-rs}{1+\theta}}\text{ in }\Omega,$

we obtain that $Tu$ satisfies

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}(Tu)&\geq&M_{2}^{-q}d(x)^{-\frac{q(2t-rs)}{1+\theta}}(Tu)^{-p},\
\ \ Tu>0&{\rm in}\ \ \Omega\\\ Tu&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right..$

Thus, $\overline{u}:=M_{2}^{\frac{q}{1+p}}Tu$ is a positive classical
supersolution to (14), because $-q+\frac{q}{1+p}=-p\frac{q}{1+p}$. By (15) and
(18) we obtain

$Tu=M_{2}^{-\frac{q}{1+p}}\overline{u}\geq
c_{1}M_{2}^{-\frac{q}{1+p}}d(x)^{s}\geq m_{1}d(x)^{s}\text{ in }\Omega.$

By inequality $m_{2}d(x)^{\frac{2t-rs}{1+\theta}}\leq v$ in $\Omega$ and the
definition of $Tu$ we conclude that

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}(Tu)&\leq&m_{2}^{-q}d(x)^{-\frac{q(2t-rs)}{1+\theta}}(Tu)^{-p},\
\ \ Tu>0&{\rm in}\ \ \Omega\\\ Tu&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right..$

Therefore, $\underline{u}:=m_{2}^{\frac{q}{1+p}}Tu$ is a positive classical
subsolution of problem (14). Hence, from (15) and (17) we have

$Tu=m_{2}^{-\frac{q}{1+p}}\overline{u}\leq
c_{2}m_{2}^{-\frac{q}{1+p}}d(x)^{s}\leq M_{1}d(x)^{s}\text{ in }\Omega.$

This way, we have proved that $Tu$ satisfies

$m_{1}d(x)^{s}\leq Tu\leq M_{1}d(x)^{s}\text{ in }\Omega.$

Similarly, using the definition of $\mathcal{A}$ and the properties of the sub
and supersolutions of problem (16) we can prove that $Tv$ satisfies

$m_{2}d(x)^{\frac{2t-rs}{1+\theta}}\leq Tv\leq
M_{2}d(x)^{\frac{2t-rs}{1+\theta}}\text{ in }\Omega.$

Then, $(Tu,Tv)\in\mathcal{A}$ for all $(u,v)\in\mathcal{A}$, that is,
$\mathcal{F}(\mathcal{A})\subseteq\mathcal{A}$.

Step 2: $\mathcal{F}$ is compact and continuous. Let $(u,v)\in\mathcal{A}$.
Then, we conclude $Tu\in C^{s}(\mathbb{R}^{n})$ and $Tv\in
C^{\frac{2t-rs}{1+\theta}}(\mathbb{R}^{n})$. Recalling that the embedding
$C^{0,s}(\overline{\Omega})\hookrightarrow
C^{0,s-\varepsilon_{1}}(\overline{\Omega})$ and
$C^{0,\frac{2t-rs}{1+\theta}}(\overline{\Omega})\hookrightarrow
C^{0,\frac{2t-rs}{1+\theta}-\varepsilon_{1}}(\overline{\Omega})$ are compact,
it follows that $\mathcal{F}$ is also compact.

Now,rest to prove that $\mathcal{F}$ is continuous. To this end, let
$(u_{n},v_{n})\subset\mathcal{A}$ be such that $u_{n}\rightarrow u$ in
$C^{s-\varepsilon_{1}}(\mathbb{R}^{n})$ and $v_{n}\rightarrow v$ in
$C^{\frac{2t-rs}{1+\theta}-\varepsilon_{1}}(\mathbb{R}^{n})$ as
$n\rightarrow\infty$. Since $\mathcal{F}$ is compact, there exists
$(U,V)\in\mathcal{A}$ such that up to a subsequence we get

$Tu_{n}\rightarrow U\text{ in }C^{s-\varepsilon_{1}}(\mathbb{R}^{n})\text{ and
}Tv_{n}\rightarrow V\text{ in
}C^{\frac{2t-rs}{1+\theta}-\varepsilon_{1}}(\mathbb{R}^{n})\text{ as
}n\rightarrow\infty.$

By Theorem 2.7 of [11], we have $(U,V)$ is a positive viscosity solution of
system (see definition in the paper [11]).

$\left\\{\begin{array}[]{llll}(-\Delta)^{s}U=v^{-q}U^{-p}&{\rm in}\ \
\Omega\\\ (-\Delta)^{t}V=u^{-r}V^{-\theta}&{\rm in}\ \ \Omega\\\ U=V=0&{\rm
in}\ \ \mathbb{R}^{n}\setminus\Omega\end{array}\right..$

From the uniqueness of positive weak solution of the problem (20), it follows
that $Tu=U$ and $Tv=V$. So,

$Tu_{n}\rightarrow Tu\text{ in }C^{s-\varepsilon_{1}}(\mathbb{R}^{n})\text{
and }Tv_{n}\rightarrow Tv\text{ in
}C^{\frac{2t-rs}{1+\theta}-\varepsilon_{1}}(\mathbb{R}^{n})\text{ as
}n\rightarrow\infty.$

So that, $\mathcal{F}$ is continuous.

Applying the Schauder’s fixed point theorem, there exists
$(u,v)\in\mathcal{A}$ such that $\mathcal{F}(u,v)=(u,v)$, that is, $Tu=u$ and
$Tv=v$. Therefore, $(u,v)$ is a positive classical solution of system (1).

The others cases will be considered similarly. But, due to the different
boundary behavior of solutions described in Proposition 2.3, the set
$\mathcal{A}$ and the constants $c_{1},c_{2}$ have to be modified accordingly.
We shall point out how these constants are chosen in order to apply the
Schauder’s fixed point theorem.

Case 2: $\frac{rs}{t}+\theta=1$ and $\alpha=\frac{qt}{s}+p<1$. By Proposition
2.3(i) and (ii) there exists $0<a<1$ and $0<c_{1}<1<c_{2}$ such that:

$\bullet$ Any positive classical subsolution $\underline{u}$ of the problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&d(x)^{-qt}u^{-p}&{\rm in}\ \
\Omega\\\ u&=&0&{\rm in}\ \ \mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

verify

$\underline{u}(x)\leq c_{2}d(x)^{s}\text{ in }\Omega.\\\ $

$\bullet$ Any positive classical supersolution $\overline{u}$ of the problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&d(x)^{-qt(t-at)}u^{-p}&{\rm
in}\ \ \Omega\\\ u&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

satisfy

$\overline{u}(x)\geq c_{1}d(x)^{s}\text{ in }\Omega.\\\ $

$\bullet$ Any positive classical subsolution $\underline{v}$ and any positive
classical supersolution $\overline{v}$ of problem (16) satisfy

$\overline{v}(x)\geq c_{1}d(x)^{t}\text{ and }\underline{v}(x)\leq
c_{2}d(x)^{t-at}\text{ in }\Omega.$

Let $\varepsilon_{1}>0$ small enough. Here $X$ stands for the Banach space

$\\{(u,v)\in C^{s-\varepsilon_{1}}(\mathbb{R}^{n})\times
C^{t-\varepsilon-\varepsilon_{1}}(\mathbb{R}^{n}):u=v=0\ {\rm in}\
\mathbb{R}^{n}\setminus\Omega\\},$

for any $\varepsilon>0$ small enough, endowed with the product norm

$\|(u,v)\|_{X}:=\|u\|_{C^{s-\varepsilon_{1}}(\mathbb{R}^{n})}+\|v\|_{C^{t-\varepsilon-\varepsilon_{1}}(\mathbb{R}^{n})}\,.$

Set

$\displaystyle\mathcal{A}:=\left\\{(u,v)\in X:m_{1}d(x)^{s}\leq u\leq
M_{1}d(x)^{s}\text{ and }m_{2}d(x)^{t}\leq v\leq M_{2}d(x)^{t-at}\text{ in
}\Omega\right\\},$

where $0<m_{i}<1<M_{i}\ (i=1,2)$ satisfy (17), (18) and

$m_{2}[diam(\Omega)]^{at}<M_{2}.$

Define the operator $\mathcal{F}$ as in the Case 1 by (20) and (21). The
inclusion $\mathcal{F}(\mathcal{A})\subseteq\mathcal{A}$ and that
$\mathcal{F}$ is continuous and compact follow as before.

Case 3: $\frac{rs}{t}+\theta<1$ and $\alpha=\frac{qt}{s}+p<1$. Let
$\varepsilon_{1}>0$ small enough. Here $X$ stands for the Banach space

$\\{(u,v)\in C^{s-\varepsilon_{1}}(\mathbb{R}^{n})\times
C^{t-\varepsilon_{1}}(\mathbb{R}^{n}):u=v=0\ {\rm in}\
\mathbb{R}^{n}\setminus\Omega\\},$

endowed with the product norm

$\|(u,v)\|_{X}:=\|u\|_{C^{s-\varepsilon_{1}}(\mathbb{R}^{n})}+\|v\|_{C^{t-\varepsilon_{1}}(\mathbb{R}^{n})}\,.$

In the same manner we define

$\displaystyle\mathcal{A}:=\left\\{(u,v)\in X:m_{1}d(x)^{s}\leq u\leq
M_{1}d(x)^{s}\text{ and }m_{2}d(x)^{t}\leq v\leq M_{2}d(x)^{t}\text{ in
}\Omega\right\\},$

where $0<m_{i}<1<M_{i}\ (i=1,2)$ satisfy (17), (18) for suitable constants
$c_{1}$ and $c_{2}$.

Case 4: $\frac{rs}{t}+\theta<1$ and $\alpha=\frac{qt}{s}+p=1$. The approach is
the same as in Case 2 above.

Let $\varepsilon_{1}>0$ small enough. Here $X$ stands for the Banach space

$\\{(u,v)\in C^{s-\varepsilon-\varepsilon_{1}}(\mathbb{R}^{n})\times
C^{t-\varepsilon_{1}}(\mathbb{R}^{n}):u=v=0\ {\rm in}\
\mathbb{R}^{n}\setminus\Omega\\},$

for any $\varepsilon>0$ small enough, endowed with the product norm

$\|(u,v)\|_{X}:=\|u\|_{C^{s-\varepsilon-\varepsilon_{1}}(\mathbb{R}^{n})}+\|v\|_{C^{t-\varepsilon_{1}}(\mathbb{R}^{n})}\,.$

Set

$\displaystyle\mathcal{A}:=\left\\{(u,v)\in X:m_{1}d(x)^{s}\leq u\leq
M_{1}d(x)^{s-as}\text{ and }m_{2}d(x)^{t}\leq v\leq M_{2}d(x)^{t}\text{ in
}\Omega\right\\},$

for some $0<a<1$, where $0<m_{i}<1<M_{i}\ (i=1,2)$ satisfy (17), (18) and

$m_{1}[diam(\Omega)]^{as}<M_{1}.$

Case 5: $\frac{rs}{t}+\theta>1$ and $\alpha=\frac{qt}{s}+p=1$. Let $0<a<1$ be
fixed such that $\frac{ars}{t}+\theta>1$. Then,

$\frac{q(2t-rs)}{s(1+\theta)}+p<1\text{ and
}\frac{q(2t-rsa)}{s(1+\theta)}+p<1.$

So, by Proposition 2.3(i), (iii), there exist $0<c_{1}<1<c_{2}$ such that:

$\bullet$ Any positive classical subsolution $\underline{u}$ of the problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&d(x)^{-\frac{q(2t-rsa)}{1+\theta}}u^{-p}&{\rm
in}\ \ \Omega\\\ u&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

verify

$\underline{u}(x)\leq c_{2}d(x)^{sa}\text{ in }\Omega.\\\ $

$\bullet$ Any positive classical supersolution $\overline{u}$ of the problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&d(x)^{-\frac{q(2t-rs)}{1+\theta}}u^{-p}&{\rm
in}\ \ \Omega\\\ u&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

satisfy

$\overline{u}(x)\geq c_{1}d(x)^{s}\text{ in }\Omega.\\\ $

$\bullet$ Any positive classical subsolution $\underline{v}$ of problem (16)
satisfies

$\underline{v}(x)\leq c_{2}d(x)^{\frac{2t-rs}{1+\theta}}\text{ in }\Omega.$

$\bullet$ Any positive classical supersolution $\overline{v}$ of problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{t}v&=&d(x)^{-ars}v^{-\theta}&{\rm
in}\ \ \Omega\\\ v&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

satisfies

$\overline{v}(x)\geq c_{1}d(x)^{\frac{2t-rsa}{1+\theta}}\text{ in }\Omega.$

Let $\varepsilon_{1}>0$ small enough. Here $X$ stands for the Banach space

$\\{(u,v)\in C^{s-\varepsilon-\varepsilon_{1}}(\mathbb{R}^{n})\times
C^{\frac{2t-rs}{1+\theta}-\varepsilon_{1}}(\mathbb{R}^{n}):u=v=0\ {\rm in}\
\mathbb{R}^{n}\setminus\Omega\\},$

for any $\varepsilon>0$ small enough, endowed with the product norm

$\|(u,v)\|_{X}:=\|u\|_{C^{s-\varepsilon-\varepsilon_{1}}(\mathbb{R}^{n})}+\|v\|_{C^{\frac{2t-rs}{1+\theta}-\varepsilon_{1}}(\mathbb{R}^{n})}\,.$

Set

$\displaystyle\mathcal{A}:=\left\\{(u,v)\in X:m_{1}d(x)^{s}\leq u\leq
M_{1}d(x)^{sa}\text{ and }m_{2}d(x)^{\frac{2t-rsa}{1+\theta}}\leq v\leq
M_{2}d(x)^{\frac{2t-rs}{1+\theta}}\text{ in }\Omega\right\\},$

where $0<m_{i}<1<M_{i}\ (i=1,2)$ satisfy (17), (18) in which the constants
$c_{1},c_{2}$ are those given above and

$m_{1}[diam(\Omega)]^{s-as}<M_{1}\text{ and
}m_{2}[diam(\Omega)]^{\frac{r(s-as)}{1+\theta}}<M_{2}.$

Case 6: $\frac{rs}{t}+\theta=1$ and $\alpha=\frac{qt}{s}+p=1$. We proceed in
the same manner as above by considering $X$ stands for the Banach space

$\\{(u,v)\in C^{s-\varepsilon-\varepsilon_{1}}(\mathbb{R}^{n})\times
C^{t-\varepsilon-\varepsilon_{1}}(\mathbb{R}^{n}):u=v=0\ {\rm in}\
\mathbb{R}^{n}\setminus\Omega\\},$

for any $\varepsilon,\varepsilon_{1}>0$ small enough, endowed with the product
norm

$\|(u,v)\|_{X}:=\|u\|_{C^{s-\varepsilon-\varepsilon_{1}}(\mathbb{R}^{n})}+\|v\|_{C^{t-\varepsilon-\varepsilon_{1}}(\mathbb{R}^{n})}\,.$

Set

$\displaystyle\mathcal{A}:=\left\\{(u,v)\in X:m_{1}d(x)^{s}\leq u\leq
M_{1}d(x)^{s-a_{1}s}\text{ and }m_{2}d(x)^{t}\leq v\leq
M_{2}d(x)^{t-a_{2}t}\text{ in }\Omega\right\\},$

where $0<a_{1},a_{2}<1$ are fixed constants and $0<m_{i}<1<M_{i}\ (i=1,2)$
satisfy (17), (18) for suitable constants $c_{1},c_{2}>0$ and

$m_{1}[diam(\Omega)]^{a_{1}s}<M_{1}\text{ and
}m_{2}[diam(\Omega)]^{a_{2}t}<M_{2}.$

(iii) Let

$a=\frac{2st\left(\frac{1+\theta}{t}-\frac{q}{s}\right)}{(1+p)(1+\theta)-qr}\text{
and }b=\frac{2st\left(\frac{1+p}{s}-\frac{r}{t}\right)}{(1+p)(1+\theta)-qr}.$

Then

$a=\frac{2s-bq}{1+p}\text{ and }b=\frac{2t-ra}{1+\theta}.$ (22)

From hypothesis, we have $a<s$ and $b<t$. Then, $bq<qt<2s$ and $ar<rs<2t$.
Now, since $\frac{bq}{s}+p>1$ and $\frac{ar}{t}+\theta>1$, from Proposition
2.3(iii) and (22) above we can find $0<c_{1}<1<c_{2}$ such that:

$\bullet$ Any positive classical subsolution $\underline{u}$ and any positive
classical supersolution $\overline{u}$ of the problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u&=&d(x)^{-bq}u^{-p}&{\rm in}\ \
\Omega\\\ u&=&0&{\rm in}\ \ \mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

satisfy

$\overline{u}(x)\geq c_{1}d(x)^{a}\text{ and }\underline{u}(x)\leq
c_{2}d(x)^{a}\text{ in }\Omega.\\\ $

$\bullet$ Any positive classical subsolution $\underline{v}$ and any positive
classical supersolution $\overline{v}$ of the problem

$\left\\{\begin{array}[]{rrll}(-\Delta)^{t}v&=&d(x)^{-ar}v^{-\theta}&{\rm in}\
\ \Omega\\\ v&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

verify

$\overline{v}(x)\geq c_{1}d(x)^{b}\text{ and }\underline{v}(x)\leq
c_{2}d(x)^{b}\text{ in }\Omega.$

As before, let $\varepsilon_{1}>0$ small enough and define $X$ to be the
Banach space

$\\{(u,v)\in C^{a-\varepsilon_{1}}(\mathbb{R}^{n})\times
C^{b-\varepsilon_{1}}(\mathbb{R}^{n}):u=v=0\ {\rm in}\
\mathbb{R}^{n}\setminus\Omega\\}$

endowed with the product norm

$\|(u,v)\|_{X}:=\|u\|_{C^{a-\varepsilon_{1}}(\mathbb{R}^{n})}+\|v\|_{C^{b-\varepsilon_{1}}(\mathbb{R}^{n})}\,.$

Set

$\displaystyle\mathcal{A}:=\left\\{(u,v)\in X:m_{1}d(x)^{a}\leq u\leq
M_{1}d(x)^{a}\text{ and }m_{2}d(x)^{b}\leq v\leq M_{2}d(x)^{b}\text{ in
}\Omega\right\\},$

where $0<m_{1}<1<M_{1}$ and $0<m_{2}<1<M_{2}$ satisfy (17) and (18). This
completes the proof of Theorem 1.2.

##  5 Proof of Theorem 1.3

We shall prove only (i); the case (ii) follows similarly.

Let $(u_{1},v_{1})$ and $(u_{2},v_{2})$ be two positive classical solutions of
system (1). Note that if $\frac{qt}{s}+p<1$, then by Theorem 1.1, we deduce
$rs<2t$. By Proposition 2.2 there exists $c_{1}>0$ such that

$u_{i}\geq c_{1}d(x)^{s}\text{ and }v_{i}\geq c_{1}d(x)^{t}$ (23)

in $\Omega$, $i=1,2$. Then, $u_{i}$ satisfies

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}u_{i}&=&v_{i}^{-q}u_{i}^{-p}\leq
c_{2}d(x)^{-qt}u_{i}^{-p},\ u_{i}>0&{\rm in}\ \ \Omega\\\ u_{i}&=&0&{\rm in}\
\ \mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

for some $c_{2}>0$. Since $\frac{qt}{s}+p<1$, by Proposition 2.3(i) and (23)
there exists $0<c<1$ such that

$cd(x)^{s}\leq u_{i}(x)\leq\frac{1}{c}d(x)^{s}$

in $\Omega$, $i=1,2$. Therefore there exists a constant $C>1$ such that
$Cu_{1}\geq u_{2}$ and $Cu_{2}\geq u_{1}$ in $\Omega$.

We claim that $u_{1}\geq u_{2}$ in $\Omega$. Supposing by contradiction, let

$\Gamma=\inf\\{\gamma>1:\gamma u_{1}\geq u_{2}\text{ in }\Omega\\}.$

By our assumption, we have $\Gamma>1$. From $\Gamma u_{1}\geq u_{2}$ in
$\Omega$, it follows that

$(-\Delta)^{t}v_{2}=u_{2}^{-r}v_{2}^{-\theta}\geq\Gamma^{-r}u_{1}^{-r}v_{2}^{-\theta}$

in $\Omega$. Thus $v_{1}$ is a positive classical solution and
$\Gamma^{\frac{r}{1+\theta}}v_{2}$ is a positive classical supersolution of

$\left\\{\begin{array}[]{rrll}(-\Delta)^{t}w&=&u_{1}^{-r}w^{-\theta}&{\rm in}\
\ \Omega\\\ w&=&0&{\rm in}\ \
\mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

because $-r+\frac{r}{1+\theta}=-\theta\frac{r}{1+\theta}$. From the
Proposition 2.1, we obtain $v_{1}\leq\Gamma^{\frac{r}{1+\theta}}v_{2}$ in
$\Omega$. Combining the above estimate, we get

$(-\Delta)^{s}u_{1}=v_{1}^{-q}u_{1}^{-p}\geq\Gamma^{-\frac{qr}{1+\theta}}v_{2}^{-q}u_{1}^{-p}$

in $\Omega$. Therefore $u_{2}$ is a positive classical solution and
$\Gamma^{\frac{qr}{(1+p)(1+\theta)}}u_{1}$ is a positive classical
supersolution of

$\left\\{\begin{array}[]{rrll}(-\Delta)^{s}z&=&v_{2}^{-q}z^{-p}&{\rm in}\ \
\Omega\\\ z&=&0&{\rm in}\ \ \mathbb{R}^{n}\setminus\Omega\end{array}\right.,$

because
$-\frac{qr}{1+\theta}+\frac{qr}{(1+p)(1+\theta)}=-p\frac{qr}{(1+p)(1+\theta)}$.
By Proposition 2.1, we conclude
$u_{2}\leq\Gamma^{\frac{qr}{(1+p)(1+\theta)}}u_{1}$ in $\Omega$. Since
$\Gamma>1$ and $\frac{qr}{(1+p)(1+\theta)}<1$, the above inequality
contradicts the minimality of $\Gamma$. Then, $u_{1}\geq u_{2}$ in $\Omega$.
Arguing similarly we conclude $u_{1}\leq u_{2}$ in $\Omega$, so $u_{1}\equiv
u_{2}$ which we obtain $v_{1}\equiv v_{2}$. Thus, the system has a unique
positive classical solution. This ends the proof of uniqueness.

## References

  * [1]
  * [2] N. Abatangelo - Large s-harmonic functions and boundary blow-up solutions for the fractional Laplacian, Discrete Contin. Dyn. Syst. 35 (2015), 5555-5607.
  * [3] Adimurthi, J. Giacomoni, S. Santra - Positive solutions to a fractional equation with singular nonlinearity, J. Differential Equations 265 (2018), 1191-1226.
  * [4] A. Ambrosetti, P. Rabinowitz - Dual variational methods in critical points theory and applications, J. Funct. Anal. 14 (1973), 349-381.
  * [5] B. Gidas, J. Spruck - Global and local behavior of positive solutions of nonlinear elliptic equations, Comm. Pure Appl. Math. 34 (1981), 525-598.
  * [6] M. Ghergu - Lane-Emden systems with negative exponents, J. Funct. Anal. 258 (2010), 3295-3318.
  * [7] E. Leite, M. Montenegro - A priori bounds and positive solutions for non-variational fractional elliptic systems, Differential and Integral Equations 30 (2017), 947-974.
  * [8] E.J.F. Leite, M. Montenegro - On positive viscosity solutions of fractional Lane-Emden systems, accepted to publication in TMNA.
  * [9] E.J.F. Leite, M. Montenegro - Principal curves to nonlocal Lane-Emden systems and related maximum principles, preprint, 2018.
  * [10] S. I. Pohozaev - On the eigenfunctions of the equation $\Delta u+\lambda f(u)=0$, Dokl. Akad. Nauk SSSR 165 (1965), 1408-1411.
  * [11] A. Quaas, A. Xia - Liouville type theorems for nonlinear elliptic equations and systems involving fractional Laplacian in the half space, Calc. Var. Partial Differential Equations 526 (2014), 1-19.
  * [12] X. Ros-Oton, J. Serra - The extremal solution for the fractional Laplacian, Calc. Var. Partial Differential Equations 50 (2014), 723-750.
  * [13] X. Ros-Oton, J. Serra - The Dirichlet problem for the fractional Laplacian: regularity up to the boundary, J. Math. Pures Appl. 101 (2014), 275-302.
  * [14] X. Ros-Oton, J. Serra - The Pohozaev identity for the fractional Laplacian, Arch. Ration. Mech. Anal. 213 (2014), 587-628.
  * [15] R. Servadei, E. Valdinoci - Variational methods for non-local operators of elliptic type, Discrete Contin. Dyn. Syst. 33 (2013), 2105-2137.
  * [16] R. Servadei, E. Valdinoci - Weak and viscosity solutions of the fractional Laplace equation, Publ. Mat. 58 (2014), 133-154.
  * [17] L. Silvestre - Regularity of the obstacle problem for a fractional power of the Laplace operator, Commun. Pure Appl. Math. 60 (2007), 67-112.


# Improving Optimal Power Flow Relaxations Using 3-Cycle Second-Order Cone
Constraints

Frederik Geth, James Foster F. Geth and J. Foster are with the Energy Systems
program, CSIRO Energy, Newcastle NSW, Australia (e-mail:
<EMAIL_ADDRESS>james.foster@csiro.au).

###### Abstract

This paper develops a novel second order cone relaxation of the semidefinite
programming formulation of optimal power flow, that does not imply the ‘angle
relaxation’. We build on a technique developed by Kim et al., extend it for
complex matrices, and apply it to $3\times 3$ positive semidefinite matrices
to generate novel second-order cone constraints that augment upon the well-
known $2\times 2$ principal-minor based second-order cone constraints.
Finally, we apply it to optimal power flow in meshed networks and provide
numerical illustrations.

###### Index Terms:

Optimal Power Flow, Mathematical Optimization

##  I Introduction and Background

Convex relaxation of nonlinear optimization models has been successfully
applied to optimal power flow in the context of certifying global optimality
or infeasibility [1], as well as in enabling mixed-integer nonlinear
optimization. The tightest known relaxations are based on semidefinite
programming (SDP) forms, either by construction, or generated as part of the
Lasserre hierarchy [2]. Since the scalability of SDP is limited in practice,
so there remains a need for tight but scalable relaxations. Second-order cone
(SOC) programming is considered more scalable, and mixed-integer (MI) SOC
solvers are commercially available.

It is noted that the physical interpretation of the difference between the SDP
and the SOC relaxation is the ‘angle relaxation’ [3], that is, Kirchhoff’s
voltage law in loops is _not_ represented beyond the implications of Ohm’s law
on voltage magnitude. Therefore, Kocuk et al. developed valid inequalities to
include the physics of 3 and 4-cycles in a SOC formulation [4], by iteratively
adding linear inequalities based on 3 and 4-cycle SDP subproblems.
Alternatively, Coffrin et al. intersect the conventional SDP/SOC BIM
power/lifted-voltage formulation with a convex hull power/polar-voltage
formulation [5], which therefore also improves upon the angle-related physics
representation. Hijazi et al. [6] develop polynomial expressions for $3\times
3$ minors and propose an interpretation w.r.t. Kirchhoff’s voltage law.

Alternatively, in this paper, we propose a method to generate tighter SOC
relaxations, based on a technique of Kim et al. [7], and apply it to SDP
optimization models for meshed power networks. It will be shown that the novel
SOC formulation improves upon the conventional, $2\times 2$ principal minor
based SOC relaxation (abbreviated PM SOC here) of the bus injection model
(BIM) optimal power flow (OPF) formulation [3]. Finally, we provide an
intuitive interpretation of the novel SOC constraints for 3-cycles by
comparing the expressions with those published by Hijazi et al. [6] and Kocuk
et al. [8].

##  II Kim et al. SDP to SOC Relaxation

Kim et al. [7] develop a SDP to SOC relaxation procedure in the real numbers,
whereas this section develops the generalization to complex matrices; staying
in the complex domain enables tighter relaxations. Specifically, we adapt the
derivation §3 ‘An extension of the SOC relaxation of type 1’, Lemma 3.2 [7] to
Hermitian matrices (a.k.a. complex symmetric).

We partition an Hermitian, positive semidefinite (PSD) variable matrix
${\color[rgb]{0,0,1}\mathbf{M}}$ as,

$\displaystyle{\color[rgb]{0,0,1}\mathbf{M}}=\begin{bmatrix}{\color[rgb]{0,0,0}\alpha}&{\color[rgb]{0,0,1}\mathbf{a}}^{\text{H}}\\\
{\color[rgb]{0,0,1}\mathbf{a}}&{\color[rgb]{0,0,1}\mathbf{A}}\\\
\end{bmatrix}\succeq 0.$ (1)

where superscript H is the conjugate transpose. Note that in general we could
position the scalar $\alpha$ anywhere on the diagonal through a permutation of
the rows and columns. This results in the following matrix variable sizes and
structures:

$\displaystyle{\color[rgb]{0,0,1}\mathbf{M}}\in{\color[rgb]{.25,.25,.25}\mathbb{H}^{n+1}},{\color[rgb]{0,0,1}\mathbf{A}}=\mathbf{A}^{\text{re}}+{\color[rgb]{1,0,0}j}\mathbf{A}^{\text{im}}\in{\color[rgb]{.25,.25,.25}\mathbb{H}^{n}},{\color[rgb]{.25,.25,.25}\mathbb{S}_{+}^{n}},$
$\displaystyle{\color[rgb]{0,0,1}\mathbf{a}}=\mathbf{a}^{\text{re}}+{\color[rgb]{1,0,0}j}\mathbf{a}^{\text{im}}\in{\color[rgb]{.25,.25,.25}\mathbb{C}^{n\times
1}},{\color[rgb]{0,0,0}\alpha}\in{\color[rgb]{.25,.25,.25}\mathbb{R}^{1\times
1}},{\color[rgb]{0,0,0}\alpha}\geq 0.$

We choose a complex matrix parameter ${\color[rgb]{.75,.5,.25}\mathbf{c}}$ and
derive ${\color[rgb]{.75,.5,.25}\mathbf{C}}$,

$\displaystyle{\color[rgb]{.75,.5,.25}\mathbf{c}}={\color[rgb]{1,0,0}\mathbf{c}}^{\text{re}}+{\color[rgb]{1,0,0}j}{\color[rgb]{1,0,0}\mathbf{c}}^{\text{im}}\in{\color[rgb]{.25,.25,.25}\mathbb{C}^{n\times
q}},\quad{\color[rgb]{.75,.5,.25}\mathbf{C}}={\color[rgb]{1,0,0}\mathbf{C}}^{\text{re}}+{\color[rgb]{1,0,0}j}{\color[rgb]{1,0,0}\mathbf{C}}^{\text{im}}={\color[rgb]{.75,.5,.25}\mathbf{c}}{\color[rgb]{.75,.5,.25}\mathbf{c}}^{\text{H}},$

or vice-versa we choose
${\color[rgb]{.75,.5,.25}\mathbf{C}}\in{\color[rgb]{.25,.25,.25}\mathbb{H}^{n}},{\color[rgb]{.25,.25,.25}\mathbb{S}_{+}^{n}}$
and factor it into ${\color[rgb]{.75,.5,.25}\mathbf{c}}$. We use $\bullet$ to
represent the Frobenius inner product of two matrices,

$\displaystyle{\color[rgb]{.75,.5,.25}\mathbf{C}}\bullet{\color[rgb]{0,0,1}\mathbf{A}}=\sum_{ij}([{\color[rgb]{.75,.5,.25}\mathbf{C}}]_{ij})^{*}[{\color[rgb]{0,0,1}\mathbf{A}}]_{ij},$
(2)

where $[\mathbf{X}]_{ij}$ means extracting the scalar at index $(i,j)$ in
matrix $\mathbf{X}$, and superscript $*$ is the complex conjugate. For PSD
Hermitian matrices this product is nonnegative,

$\displaystyle{\color[rgb]{.75,.5,.25}\mathbf{C}}\bullet{\color[rgb]{0,0,1}\mathbf{A}}\geq
0.$ (3)

We now illustrate Kim et al.’s main result, a valid convex quadratic
inequality for a sampled
${\color[rgb]{.75,.5,.25}\mathbf{C}}={\color[rgb]{.75,.5,.25}\mathbf{c}}{\color[rgb]{.75,.5,.25}\mathbf{c}}^{\text{H}}$,
for any Hermitian and PSD matrix ${\color[rgb]{0,0,1}\mathbf{M}}$ partitioned
with a scalar $\alpha$ on the diagonal:

$\displaystyle({\color[rgb]{.75,.5,.25}\mathbf{c}}^{\text{H}}{\color[rgb]{0,0,1}\mathbf{a}})^{\text{H}}({\color[rgb]{.75,.5,.25}\mathbf{c}}^{\text{H}}{\color[rgb]{0,0,1}\mathbf{a}})\leq{\color[rgb]{0,0,0}\alpha}\cdot({\color[rgb]{.75,.5,.25}\mathbf{C}}\bullet{\color[rgb]{0,0,1}\mathbf{A}}).$
(4)

Implementation can use identities w.r.t. the reals,

$\displaystyle{\color[rgb]{.75,.5,.25}\mathbf{C}}\bullet{\color[rgb]{0,0,1}\mathbf{A}}=\mathbf{C}^{\text{re}}\bullet\mathbf{A}^{\text{re}}+\mathbf{C}^{\text{im}}\bullet\mathbf{A}^{\text{im}},$
$\displaystyle\operatorname{Re}({\color[rgb]{.75,.5,.25}\mathbf{c}}^{\text{H}}{\color[rgb]{0,0,1}\mathbf{a}})=({\color[rgb]{1,0,0}\mathbf{c}}^{\text{re}})^{\text{T}}\mathbf{a}^{\text{re}}+({\color[rgb]{1,0,0}\mathbf{c}}^{\text{im}})^{\text{T}}\mathbf{a}^{\text{im}},$
$\displaystyle\operatorname{Im}({\color[rgb]{.75,.5,.25}\mathbf{c}}^{\text{H}}{\color[rgb]{0,0,1}\mathbf{a}})=({\color[rgb]{1,0,0}\mathbf{c}}^{\text{re}})^{\text{T}}\mathbf{a}^{\text{im}}-({\color[rgb]{1,0,0}\mathbf{c}}^{\text{im}})^{\text{T}}\mathbf{a}^{\text{re}}.$

Note that we directly derive (real-value) SOCs for complex-value PSD
constraints, _without_ having to invoke the isomorphism between the complex
numbers and $2\times 2$ matrices (in its block matrix form) on
${\color[rgb]{0,0,1}\mathbf{A}}$ or ${\color[rgb]{0,0,1}\mathbf{M}}$.

##  III Application to a $3\times 3$ PSD Voltage Matrix

We focus on a 3-bus fully connected power network in this section, and discuss
generalization later. Using the BIM OPF [3] model for a 3-bus system, with bus
voltages
${\color[rgb]{0,0,1}U_{i}}=U^{\text{re}}_{i}+{\color[rgb]{1,0,0}j}U^{\text{im}}_{i}$,
the system ‘${\color[rgb]{0,0,1}\mathbf{M}}$’ matrix in the ‘W’ variable space
corresponds to:

$\displaystyle\begin{bmatrix}U^{\text{re}}_{1}+{\color[rgb]{1,0,0}j}U^{\text{im}}_{1}\\\
U^{\text{re}}_{2}+{\color[rgb]{1,0,0}j}U^{\text{im}}_{2}\\\
U^{\text{re}}_{3}+{\color[rgb]{1,0,0}j}U^{\text{im}}_{3}\\\
\end{bmatrix}\begin{bmatrix}U^{\text{re}}_{1}+{\color[rgb]{1,0,0}j}U^{\text{im}}_{1}\\\
U^{\text{re}}_{2}+{\color[rgb]{1,0,0}j}U^{\text{im}}_{2}\\\
U^{\text{re}}_{3}+{\color[rgb]{1,0,0}j}U^{\text{im}}_{3}\\\
\end{bmatrix}^{\text{H}}$ (5)
$\displaystyle=\begin{bmatrix}W^{\text{re}}_{11}&W^{\text{re}}_{12}+{\color[rgb]{1,0,0}j}W^{\text{im}}_{12}&W^{\text{re}}_{13}+{\color[rgb]{1,0,0}j}W^{\text{im}}_{13}\\\
W^{\text{re}}_{12}-{\color[rgb]{1,0,0}j}W^{\text{im}}_{12}&W^{\text{re}}_{22}&W^{\text{re}}_{23}+{\color[rgb]{1,0,0}j}W^{\text{im}}_{23}\\\
W^{\text{re}}_{13}-{\color[rgb]{1,0,0}j}W^{\text{im}}_{13}&W^{\text{re}}_{23}-{\color[rgb]{1,0,0}j}W^{\text{im}}_{23}&W^{\text{re}}_{33}\\\
\end{bmatrix}={\color[rgb]{0,0,1}\mathbf{M}}\succeq 0$

The partitions (1) then are,

$\displaystyle{\color[rgb]{0,0,1}\mathbf{A}}\\!:\\!\begin{bmatrix}W^{\text{re}}_{22}\\!\\!&\\!\\!W^{\text{re}}_{23}\\!+\\!{\color[rgb]{1,0,0}j}W^{\text{im}}_{23}\\\
W^{\text{re}}_{23}\\!-\\!{\color[rgb]{1,0,0}j}W^{\text{im}}_{23}\\!\\!&\\!\\!W^{\text{re}}_{33}\\\
\end{bmatrix},{\color[rgb]{0,0,1}\mathbf{a}}\\!:\\!\begin{bmatrix}W^{\text{re}}_{12}\\!-\\!{\color[rgb]{1,0,0}j}W^{\text{im}}_{12}\\\
W^{\text{re}}_{13}\\!-\\!{\color[rgb]{1,0,0}j}W^{\text{im}}_{13}\end{bmatrix},\alpha\\!:\\!W^{\text{re}}_{11}.$

###  III-A $2\times 2$ Principal Minor SOC Constraints

For this partition, substituting
${\color[rgb]{.75,.5,.25}\mathbf{c}}=\begin{bmatrix}1&0\end{bmatrix}^{\text{T}}$
or $\begin{bmatrix}0&1\end{bmatrix}^{\text{T}}$ into (4) generates known
$2\times 2$ PM SOCs, i.e.

$\displaystyle(W^{\text{re}}_{12})^{2}+(W^{\text{im}}_{12})^{2}\leq
W^{\text{re}}_{11}W^{\text{re}}_{22}.$ (6)
$\displaystyle(W^{\text{re}}_{13})^{2}+(W^{\text{im}}_{13})^{2}\leq
W^{\text{re}}_{11}W^{\text{re}}_{33},$ (7a)

A partition with $\alpha=W^{\text{re}}_{22}$ or $\alpha=W^{\text{re}}_{33}$
and a one-hot ${\color[rgb]{.75,.5,.25}\mathbf{c}}$ generates the third unique
PM SOC, i.e.,

$\displaystyle(W^{\text{re}}_{23})^{2}+(W^{\text{im}}_{23})^{2}\leq
W^{\text{re}}_{22}W^{\text{re}}_{33}.$ (8a)

###  III-B Kim et al. SOC Constraints

We sample from the space of vectors
${\color[rgb]{.75,.5,.25}\mathbf{c}}\in{\color[rgb]{.25,.25,.25}\mathbb{C}^{2\times
1}}$ with the first element normalized, and parameterize through
${\color[rgb]{1,0,0}r},{\color[rgb]{1,0,0}\theta}$,

$\displaystyle{\color[rgb]{.75,.5,.25}\mathbf{c}}=\begin{bmatrix}1\\\
{\color[rgb]{1,0,0}r}e^{{\color[rgb]{1,0,0}j}{\color[rgb]{1,0,0}\theta}}\end{bmatrix}=\begin{bmatrix}1\\\
{\color[rgb]{1,0,0}r}\cos({\color[rgb]{1,0,0}\theta})\end{bmatrix}+{\color[rgb]{1,0,0}j}\begin{bmatrix}0\\\
{\color[rgb]{1,0,0}r}\sin({\color[rgb]{1,0,0}\theta})\end{bmatrix},$

Next we obtain a parameterized family of valid SOCs:

$\displaystyle(W^{\text{re}}_{12}+{\color[rgb]{1,0,0}r}\cos({\color[rgb]{1,0,0}\theta})W^{\text{re}}_{13}-{\color[rgb]{1,0,0}r}\sin({\color[rgb]{1,0,0}\theta})W^{\text{im}}_{13})^{2}$
$\displaystyle+(W^{\text{im}}_{12}+{\color[rgb]{1,0,0}r}\cos({\color[rgb]{1,0,0}\theta})W^{\text{im}}_{13}+{\color[rgb]{1,0,0}r}\sin({\color[rgb]{1,0,0}\theta})W^{\text{re}}_{13})^{2}$
$\displaystyle\\!\\!\leq
W^{\text{re}}_{11}\left(W^{\text{re}}_{22}+{\color[rgb]{1,0,0}r}^{2}W^{\text{re}}_{33}+2{\color[rgb]{1,0,0}r}\\!\left(\cos({\color[rgb]{1,0,0}\theta})W^{\text{re}}_{23}-\sin({\color[rgb]{1,0,0}\theta})W^{\text{im}}_{23}\right)\right)\\!.$
(9)

We also generate them for the other partitions, i.e.,

$\displaystyle(W^{\text{re}}_{12}+{\color[rgb]{1,0,0}r}\cos({\color[rgb]{1,0,0}\theta})W^{\text{re}}_{23}-{\color[rgb]{1,0,0}r}\sin({\color[rgb]{1,0,0}\theta})W^{\text{im}}_{23})^{2}$
$\displaystyle+(W^{\text{im}}_{12}-{\color[rgb]{1,0,0}r}\cos({\color[rgb]{1,0,0}\theta})W^{\text{im}}_{23}-{\color[rgb]{1,0,0}r}\sin({\color[rgb]{1,0,0}\theta})W^{\text{re}}_{23})^{2}$
$\displaystyle\\!\\!\leq
W^{\text{re}}_{22}\left(W^{\text{re}}_{11}+{\color[rgb]{1,0,0}r}^{2}W^{\text{re}}_{33}+2{\color[rgb]{1,0,0}r}\\!\left(\cos({\color[rgb]{1,0,0}\theta})W^{\text{re}}_{13}-\sin({\color[rgb]{1,0,0}\theta})W^{\text{im}}_{13}\right)\right)\\!.$
(10)
$\displaystyle(W^{\text{re}}_{13}+{\color[rgb]{1,0,0}r}\cos({\color[rgb]{1,0,0}\theta})W^{\text{re}}_{23}+{\color[rgb]{1,0,0}r}\sin({\color[rgb]{1,0,0}\theta})W^{\text{im}}_{23})^{2}$
$\displaystyle+(W^{\text{im}}_{13}+{\color[rgb]{1,0,0}r}\cos({\color[rgb]{1,0,0}\theta})W^{\text{im}}_{23}-{\color[rgb]{1,0,0}r}\sin({\color[rgb]{1,0,0}\theta})W^{\text{re}}_{23})^{2}$
$\displaystyle\\!\\!\\!\\!\leq
W^{\text{re}}_{33}\\!\left(W^{\text{re}}_{11}+{\color[rgb]{1,0,0}r}^{2}W^{\text{re}}_{22}+2{\color[rgb]{1,0,0}r}\\!\left(\cos({\color[rgb]{1,0,0}\theta})W^{\text{re}}_{12}-\sin({\color[rgb]{1,0,0}\theta})W^{\text{im}}_{12}\right)\right)\\!.$
(11)

Finally, we still need to state the linear inequalities based on the Frobenius
inner product (3), i.e.,

$\displaystyle
W^{\text{re}}_{22}+{\color[rgb]{1,0,0}r^{2}}W^{\text{re}}_{33}+2{\color[rgb]{1,0,0}r}\left(\cos({\color[rgb]{1,0,0}\theta})W^{\text{re}}_{23}-\sin({\color[rgb]{1,0,0}\theta})W^{\text{im}}_{23}\right)\geq
0,$ (12a) $\displaystyle
W^{\text{re}}_{11}+{\color[rgb]{1,0,0}r^{2}}W^{\text{re}}_{33}+2{\color[rgb]{1,0,0}r}\left(\cos({\color[rgb]{1,0,0}\theta})W^{\text{re}}_{13}-\sin({\color[rgb]{1,0,0}\theta})W^{\text{im}}_{13}\right)\geq
0,$ (13a) $\displaystyle
W^{\text{re}}_{11}+{\color[rgb]{1,0,0}r^{2}}W^{\text{re}}_{22}+2{\color[rgb]{1,0,0}r}\left(\cos({\color[rgb]{1,0,0}\theta})W^{\text{re}}_{12}-\sin({\color[rgb]{1,0,0}\theta})W^{\text{im}}_{12}\right)\geq
0.$ (14a)

Inequalities (9)-(11) are families of SOC constraints. One re-obtains a PM SOC
for ${\color[rgb]{1,0,0}r}=0$, and similarly if one scales the second element
of ${\color[rgb]{.75,.5,.25}\mathbf{c}}$ instead.

##  IV Interpretation with Kirchhoff’s Voltage Law

###  IV-A Kirchhoff’s Voltage Law for 3-Cycle in Lifted Voltage

Hijazi et al. [6] show that Kirchhoff’s voltage law (KVL) for a 3-cycle in the
lifted voltage variable space is,

$\displaystyle(W^{\text{re}}_{13}+{\color[rgb]{1,0,0}j}W^{\text{im}}_{13})^{*}(W^{\text{re}}_{23}+{\color[rgb]{1,0,0}j}W^{\text{im}}_{23})={W^{\text{re}}_{33}}(W^{\text{re}}_{12}+{\color[rgb]{1,0,0}j}W^{\text{im}}_{12})^{*}.$
(15)

We note that this corresponds to a minor that involves only a single diagonal
element, i.e. ‘type 2’ in [8]. In the reals, this becomes,

$\displaystyle
W^{\text{re}}_{13}W^{\text{re}}_{23}+W^{\text{im}}_{13}W^{\text{im}}_{23}=W^{\text{re}}_{33}W^{\text{re}}_{12},$
(16a) $\displaystyle
W^{\text{re}}_{13}W^{\text{im}}_{23}-W^{\text{im}}_{13}W^{\text{re}}_{23}=-W^{\text{re}}_{33}W^{\text{im}}_{12}.$
(17a)

The real-value permuted 3-cycle KVL expressions are,

$\displaystyle
W^{\text{re}}_{12}W^{\text{re}}_{13}-W^{\text{im}}_{12}W^{\text{im}}_{13}=W^{\text{re}}_{11}W^{\text{im}}_{23},$
(18a) $\displaystyle
W^{\text{re}}_{12}W^{\text{im}}_{13}+W^{\text{im}}_{12}W^{\text{re}}_{13}=W^{\text{re}}_{11}W^{\text{im}}_{23},$
(19a) $\displaystyle
W^{\text{re}}_{12}W^{\text{re}}_{23}+W^{\text{im}}_{12}W^{\text{im}}_{23}=W^{\text{re}}_{22}W^{\text{re}}_{13},$
(20a) $\displaystyle
W^{\text{re}}_{12}W^{\text{im}}_{23}+W^{\text{im}}_{12}W^{\text{re}}_{23}=-W^{\text{re}}_{22}W^{\text{im}}_{13}.$
(21a)

Figure 1: Contour plot for ‘case3_lmbd’

Figure 2: Contour plot for ‘case5_pjm’

Figure 3: Contour plot for ‘case14_ieee’

###  IV-B Kim et al. SOCs with Solution in Equality

If $\operatorname{rank}({\color[rgb]{0,0,1}\mathbf{M}})=1$, PM SOCs (7a)-(8a)
are active,

$\displaystyle(W^{\text{re}}_{13})^{2}+(W^{\text{im}}_{13})^{2}=W^{\text{re}}_{33}W^{\text{re}}_{11},$
(22a)
$\displaystyle(W^{\text{re}}_{23})^{2}+(W^{\text{im}}_{23})^{2}=W^{\text{re}}_{33}W^{\text{re}}_{22}.$
(23a)

Under the same conditions, substituting
${\color[rgb]{1,0,0}r_{1}}={\color[rgb]{1,0,0}r_{2}}=1,$ and
${\color[rgb]{1,0,0}\theta_{1}}=0,{\color[rgb]{1,0,0}\theta_{2}}=\pi/2$, in
(9) we obtain,

$\displaystyle(W^{\text{re}}_{13}+W^{\text{re}}_{23})^{2}+(W^{\text{im}}_{13}+W^{\text{im}}_{23})^{2}$
$\displaystyle=W^{\text{re}}_{33}(W^{\text{re}}_{11}+W^{\text{re}}_{22}+2W^{\text{re}}_{12}),$
(24a)
$\displaystyle(W^{\text{re}}_{13}+W^{\text{im}}_{23})^{2}+(W^{\text{im}}_{13}-W^{\text{re}}_{23})^{2}$
$\displaystyle=W^{\text{re}}_{33}(W^{\text{re}}_{11}+W^{\text{re}}_{22}-2W^{\text{im}}_{12}).$
(25a)

Therefore the following identities hold when
$\operatorname{rank}({\color[rgb]{0,0,1}\mathbf{M}})=1$,

$\displaystyle\eqref{eq_pm_ccaa}+\eqref{eq_pm_ccbb}+2\cdot\eqref{eq_cycle_c_real}\equiv\eqref{eq_kk_re_cc},$
$\displaystyle\eqref{eq_pm_ccaa}+\eqref{eq_pm_ccbb}+2\cdot\eqref{eq_cycle_c_imag}\equiv\eqref{eq_kk_im_cc},$

as well as the equivalents for (18a)-(21a) (i.e. a total of 6 SOCs). This
establishes that the Kim et al. SOCs include a stronger representation of KVL
than just the PM SOCs. Note that this _completing the square_ approach
actually embeds nonconvex equalities (16a)-(21a) into a convex quadratic form.

Note that by themselves, (24a)-(25a) do not prove (16a)-(17a) hold;
(22a)-(23a) are needed. Therefore observing that (9)-(11) are active in the
solution is insufficient to prove that
$\operatorname{rank}({\color[rgb]{0,0,1}\mathbf{M}})=1$. The PM SOCs
(22a)-(23a) must also be active in the solution to give that confidence.

##  V Numerical Illustration

We extend the BIM implementation in PowerModels [9] and use 3-cycle PGLib test
cases [10] with two ${\color[rgb]{1,0,0}\theta}$-parameterized SOCs per
3-cycle permutation. We use Ipopt [11] as the solver for NLP ‘AC’ OPF; Hypatia
[12] for the BIM SDP and the Kim et al. SOC forms, and again Ipopt for the PM
SOC. Fig. 3 \- Fig. 3 show sweeps for two SOC families
${\color[rgb]{1,0,0}\theta_{1}},{\color[rgb]{1,0,0}\theta_{2}}$, w.r.t.
objective value, for the 3, 5 and 14 bus test cases. The optimal
${\color[rgb]{1,0,0}\theta}$ for the 3-bus case is around $3.7$, with an
objective value of $5790$. For ‘case5_pjm’, we want two distinct angles, i.e.
${\color[rgb]{1,0,0}\theta_{1}}\approx
1.6,{\color[rgb]{1,0,0}\theta_{2}}\approx 4.9$ with an objective of $16181$;
for ‘case14_ieee’, one suffices, i.e. any
${\color[rgb]{1,0,0}\theta}\in[\pi;2\pi]$ gives objective $2178$. It is
interesting that, to obtain the tightest relaxation, one case requires two
fine-tuned angles, another requires one fine-tuned angle, and the third only
requires one angle in a broader range. A condition for recognizing that a
single angle is sufficient is that the maximum is on the diagonal of the plot.
Table I summarizes numerical results in comparison with the NLP, SDP and PM
SOC relaxations.

TABLE I: Objective values; Kim et al. SOCs for ${\color[rgb]{1,0,0}\theta_{1}}=0,{\color[rgb]{1,0,0}\theta_{2}}=3\pi/2$. | AC NLP | BIM SDP | Kim+PM SOC | PM SOC
---|---|---|---|---
| $/h | % gap | % gap | % gap
case3_lmbd | 5812.64 | 0.38 | 0.54 | 1.32
case5_pjm | 17551.89 | 5.22 | 14.47 | 14.54
case14_ieee | 2178.08 | 0.00 | 0.00 | 0.11

We note that in general, for the convergence of the interior point method to
be strong, you want to avoid having a large amount of inequalities active in
the solution. In general we expect this to be the case as many OPF problems
are naturally low-rank [13]. Therefore, it may be beneficial to use the Kim et
al. SOCs algorithmically as nonlinear cuts.

##  VI Conclusions

We develop a SOC relaxation does not entail the ‘angle relaxation’. For a
$3\times 3$ Hermitian matrix, the Kim et al. SOCs sampled in two directions
each ($2\times 3$ unique ones), when combined with the $2\times 2$ principal
minor-derived SOCs (3 unique ones), clearly offers an tighter SOC relaxation.
Combining with results by Hijazi et al. [6], we know now that when all 9 SOCs
are satisfied in equality, it furthermore certifies that the original matrix
is PSD and rank-1, and that the power flow satisfies Kirchhoff’s voltage law
in a 3-cycle.

Nevertheless, these additional 6 SOCs may not be sufficient to obtain a
relaxation that is equivalent to the SDP relaxation, as there are actually
families of parameteric SOCs from which one needs to sample. The results do
not suggest an obvious sampling strategy. Profiling on benchmark libraries
[10] can aid in strategy development.

Real transmissions networks have loops in the topology that go beyond
3-cycles. Nevertheless, using chordal completion similar to [4], [14], we can
re-cast the network graph in terms of intersected 3-cycles, and can apply the
technique as such. It would be interesting to combine the Kim et al. SOCs with
the QC formulation [5], to obtain a more scalable form than QC+SDP, and more
accurate than QC+PM SOC. In future work, we plan to explore generating
relaxations to $4\times 4$ and higher-dimensional PSD variable matrices
directly. A prime use case is _unbalanced_ SDP BIM OPF, as it relies on a
large amount of $6\times 6$ complex PSD matrices, and solvers struggle with
scalability and numerical stability [15].

## References

  * [1] S. Gopinath, H. L. Hijazi, T. Weisser, H. Nagarajan, M. Yetkin, K. Sundar, and R. W. Bent, “Proving global optimality of ACOPF solutions,” _Electric Power Syst. Res._ , vol. 189, no. October 2019, p. 106688, 2020\.
  * [2] C. Josz and D. K. Molzahn, “Lasserre hierarchy for large scale polynomial optimization in real and complex variables,” _SIAM Journal on Optimization_ , vol. 28, no. 2, pp. 1017–1048, 2018.
  * [3] S. H. Low, “Convex relaxation of optimal power flow - part I: formulations and equivalence,” _IEEE Trans. Control Netw. Syst._ , vol. 1, no. 1, pp. 15–27, mar 2014.
  * [4] B. Kocuk, S. Dey, and S. Andy, “Strong SOCP relaxations of optimal power flow,” _Operations Res._ , vol. 64, no. 6, pp. 1177–1196, apr 2016.
  * [5] C. Coffrin, H. Hijazi, and P. Van Hentenryck, “The QC relaxation: a theoretical and computational study on optimal power flow,” _IEEE Trans. Power Syst._ , vol. 31, no. 4, pp. 3008–3018, 2016.
  * [6] H. Hijazi, C. Coffrin, and P. Van Hentenryck, “Polynomial SDP cuts for optimal power flow,” in _PSCC_ , Genoa, Italy, 2016.
  * [7] S. Kim, M. Kojima, and M. Yamashita, “Second order cone programming relaxation of a positive semidefinite constraint,” _Optimization Methods Software_ , vol. 18, no. 5, pp. 535–541, 2003.
  * [8] B. Kocuk, S. S. Dey, and X. A. Sun, “Matrix minor reformulation and SOCP-based spatial branch-and-cut method for the AC optimal power flow problem,” 2017. [Online]. Available: https://arxiv.org/pdf/1703.03050.pdf
  * [9] C. Coffrin, R. Bent, K. Sundar, Y. Ng, and M. Lubin, “PowerModels.jl: an open-source framework for exploring power flow formulations,” in _Power Syst. Comp. Conf._ , vol. 20, Dublin, Ireland, 2018, p. 8.
  * [10] S. Babaeinejadsarookolaee, A. Birchfield, R. D. Christie, C. Coffrin, C. DeMarco, R. Diao, M. Ferris, S. Fliscounakis, S. Greene, R. Huang, C. Josz, R. Korab, B. Lesieutre, J. Maeght, D. K. Molzahn, T. J. Overbye, P. Panciatici, B. Park, J. Snodgrass, and R. Zimmerman, “The power grid library for benchmarking AC optimal power flow algorithms,” _[math.OC]_ , pp. 1–17, 2019.
  * [11] A. Wächter and L. T. Biegler, “On the implementation of primal-dual interior point filter line search algorithm for large-scale nonlinear programming,” _Math. Prog._ , vol. 106, no. 1, pp. 25–57, 2006.
  * [12] C. Coey, L. Kapelevich, and J. P. Vielma, “Towards practical generic conic optimization,” pp. 1–25, 2020. [Online]. Available: http://arxiv.org/abs/2005.01136
  * [13] J. Lavaei and S. H. Low, “Zero duality gap in optimal power flow problem,” _IEEE Trans. Power Syst._ , vol. 27, no. 1, pp. 92–107, 2012.
  * [14] R. A. Jabr, “Exploiting sparsity in SDP relaxations to the OPF problem,” _IEEE Trans. Power Syst._ , vol. 27, no. 2, pp. 1138–1139, 2012.
  * [15] L. Gan and S. H. Low, “Convex relaxations and linear approximation for optimal power flow in multiphase radial networks,” in _Power Syst. Comp. Conf._ , Wroclaw, Poland, 2014, pp. 1–9.