\section{Introduction to the Hardy-Littlewood-Sobolev type systems}
The well known Hardy-Littlewood-Sobolev inequality (HLS) states:
\begin{align}
    \int_{R^n}\int_{R^n} \displaystyle{\frac{f(x)g(y)}{|x-y|^{n-\gamma}}} dx dy
    \leq
    C(n,s,\gamma) ||f||_r||g||_s
 \end{align}
where $0 <\gamma < n$, $1< s, r  < \infty$, $\frac{1}{r} +
\frac{1}{s} + \frac{n-\gamma}{n}= 2$, $f \in L^r(R^n)$ and $g \in
L^s(R^n)$.

Define an operator $T$ such that $ Tg(x):=\displaystyle{\int}_{R^n} \displaystyle{\frac{g(y)}{|x-y|^{n-\gamma}} }dy$, $\gamma \in (0,n)$, then the HLS inequality becomes:
\begin{align}
	 ||Tg||_{\frac{ns}{n-s\gamma}} \leq C(n,s,\gamma) ||g||_s, \text{ or }
	 ||Tg||_p \leq C(n,s,\gamma) ||g||_{\frac{np}{n+\gamma p}},
\end{align}
where $\frac{n}{n-\gamma}<p< \infty$, and $1<s< n/{\gamma}$.

The best constant $C = C(n,s,\gamma)$ is the maximum of:
\begin{equation}\label{optimizingProblem}
    J(f,g) =  \int_{R^n}\int_{R^n}\frac{f(x)g(y)}{|x-y|^{n-\gamma}} dx dy
\end{equation}
with constraints $ \|f\|_r = \|g\|_s =1$.
The above optimizing problem leads us to a system of integral equations on $f$ and $g$.
Let $u = c_1 f^{r-1}$, $v = c_2 g^{s-1}$,
 $p=\frac{1}{r-1}$,
 $q=\frac{1}{s-1}$, and  choose suitable constants $c_1$ and
 $c_2$, we arrive at the
following system of Euler-Lagrange equations for the HLS
inequality :

\begin{equation} \left \{
 \begin{array}{l}
      u(x) = \int_{R^{n}} \frac{ v^q(y)}{|x-y|^{n-\gamma}} dy\\
      v(x) = \int_{R^{n}} \frac{ u^p(y)}{|x-y|^{n-\gamma}} dy
\end{array}
\right. \label{sys}
\end{equation}
with $ u, v >0$, $\;u \in
L^{p+1}$, $\;v \in L^{q+1}$, $0 < p < \infty$, $\;0 < q < \infty$,
$\;\frac{1}{p+1} + \frac{1}{q+1} =  \frac{n-\gamma}{n}$.

For the same $p,q,n,\gamma$, the solutions of  (\ref{sys}) also are solutions of (the reverse is nontrivial and might be false \cite{CLO05}):
\begin{equation} \left \{
 \begin{array}{l}
      (-\Delta)^{\gamma/2} u = v^q , \;\; u>0, \mbox{ in } R^n,\\
      (-\Delta)^{\gamma/2} v = u^p, \;\; v>0, \mbox{ in } R^n.
\end{array}
\right. \label{pde}
\end{equation}
For $0 < p, q < \infty$, we call \eqref{sys} and \eqref{pde} \textbf{Hardy-Littlewood-Sobolev (HLS) type systems}. This class of systems include some well-known special cases, for instance,
if $ p=q=\frac{n+\gamma}{n-\gamma}$, and $u(x)=v(x)$, the above system
becomes:

\begin{equation}
(-\Delta)^{\gamma/2} u = u^{(n+\gamma)/(n-\gamma)},\;\; u>0, \mbox{ in } R^n.
\label{pde10}
\end{equation}
In particular, when $n \geq 3 $, and $\gamma = 2,$
\begin{equation}
 -\Delta u =  u^{ (n+2)/(n-2) }, \;\; u>0, \mbox{
in } R^n. \label{pde2}
\end{equation}
Equation (\ref{pde10}) is equivalent to:
\begin{equation}
    u(x) = \int_{R^{n}} \frac{u(y)^{\frac{n+\gamma}{n-\gamma} }}{|x - y|^{n - \gamma}}
      dy,\;\; u >0 \mbox{ in } R^n,
      \label{sys10}
\end{equation}
for which we refer readers to Chen, Li and Ou \cite{CLO06}, where \eqref{pde10} is defined for non-integer $\gamma$.

The HLS type systems can be categorized into three cases, namely, critical case
$\frac{1}{p+1} + \frac{1}{q+1} =  \frac{n-\gamma}{n}$, subcritical case $\frac{1}{p+1} +\frac{1}{q+1} > \frac{{n-\gamma}}{n}$, supercritical case $\frac{1}{p+1} +\frac{1}{q+1} < \frac{{n-\gamma}}{n}$. The HLS type systems behave very differently in each case, hence different methods are needed to study them.

 \textbf{The critical} case of the HLS type system, also known as the HLS system, was studied in an elegant paper of E. Lieb \cite{Lieb83}. The existence of solution is proved (see \cite{Lieb83} and details in Section 3), so the remaining question is the classification/uniqueness of positive solutions.
We conjecture that all positive solutions of (\ref{sys}) or (\ref{pde}) are given by a single
family of solutions via translation and scaling.
In other words, up to some translations and scalings,
the positive solutions of (\ref{sys}) or (\ref{pde}) are unique. This kind of uniqueness is called essential uniqueness. This problem is solved in some special cases but is still open in general.

In \textbf{the supercritical case}\,
the existence of solution is also established (see Section 3).
The existence proof is based on a relatively new method which combines shooting method with degree theory. We give an outline of the method to prove the existence for an even more general system in Section 4.1, which contains supercrtical and critical HLS type systems as special cases.

We are also interested in asymptotic analysis of the solutions of supercritical HLS type systems. For example, what are the asymptotic expansions of the radial solutions?
Are all radial solutions scaling related?
These questions are addressed in Section 4.2.

Last, we consider \textbf{the subcritical} HLS type systems, in particular, the Lane-Emden system.
The so-called Lane-Emden conjecture states that,
for $0<p,q<\infty$, $\frac{1}{p+1} +\frac{1}{q+1} > \frac{n-2}{n}$,
\begin{equation}
 \left \{
 \begin{array}{l}
      -\Delta u(x) = v^q(x) , \;\; u\geq0, \mbox{ in } R^n,\\
      -\Delta v(x) = u^p(x) , \;\; v\geq0, \mbox{ in } R^n,
\end{array}
\right. \label{LaneEmden}
\end{equation}
has $u=0$ and $v=0$ as the unique locally bounded solution.

The conjecture naturally generalize to the systems \eqref{sys} or \eqref{pde} in the subcritical cases with an additional condition,
$pq>1$. Notice that $pq>1$ is a necessary condition for this conjecture to hold in high order HLS type systems. For example if $p=q=1$ and $\gamma=4$ we have solution $u=v=e^{w\cdot x}$ to \eqref{pde} for $w\in\mathbb{R}^n$ with $|w|=1$.

This paper is organized as following. In Section 2 we provide some useful estimates. In Section 3, we discuss about the classification of solution to critical HLS type system. In Section 4, existence and asymptotic analysis of solutions to supercritical cases are investigated. In Section 5, we focus on subcritical cases, in particular, the Lane-Emden system. Throughout this paper, we assume that$0<p,q<\infty$ and $pq>1$, and all positive solution are assumed locally bounded, unless specified otherwise.

\section{Some basic estimates}
For $pq>1$, denote the scaling component of system \eqref{LaneEmden} by
	\begin{align}
	\alpha = \dfrac{2(p+1)}{pq-1}, \ \ \beta = \dfrac{2(q+1)}{pq-1}.
	\end{align}

	\begin{remark} For $pq>1$, the critical hyperbola has a new form in terms of $\alpha$ and $\beta$,
		\begin{align*}
			\frac{1}{p+1} +\frac{1}{q+1} =
			\frac{n-2}{n} \Leftrightarrow \alpha +\beta = n-2.
		\end{align*}
		Hence, supercritical condition $\frac{1}{p+1} +\frac{1}{q+1} <
					\frac{n-2}{n}\Leftrightarrow \alpha +\beta < n-2$, and subcritical condition $\frac{1}{p+1} +\frac{1}{q+1} >
										\frac{n-2}{n}\Leftrightarrow \alpha +\beta > n-2$.
	\end{remark}
Here we present some useful estimates including comparison principle and energy estimates for \eqref{LaneEmden}. These estimates are valid for all three cases, i.e., critical, supercritical and subcritical. They are useful in many aspects, such as in asymptotic analysis for solutions to critical and supercritical cases, and to prove Liouville type theorem in subcritical case (Actually, to prove Lane-Emden conjecture a proper energy estimate is the key, which is detailed in Section 5).
These estimates are:
\begin{lemma}\label{energyEstimates}
Let $p,q>0$ with $pq>1$. For any positive (locally bounded) solution $(u,v)$ of \eqref{LaneEmden}
\begin{align}\label{uvEstimate1}
	\int_{B_R} u \leq C R^{n-\alpha} \ \text{and } \int_{B_R} v \leq C R^{n-\beta},
\end{align}
\begin{align}\label{uvEsitmatePQ}
	\int_{B_R} u^q \leq C R^{n-q\alpha} \ \text{and } \int_{B_R} v^p \leq C R^{n-p\beta},
\end{align}
\end{lemma}
and by Maximum-principle we get,
\begin{lemma}[Comparison principle]\label{comparisonPrinciple}
	Let $p\geq q >0$ with $pq>1$. Let $(u,v)$ be a positive bounded solution of \eqref{LaneEmden}. Then we have the following comparison principle,
	\begin{align*}
		\frac{v^{p+1}}{p+1}\leq \frac{u^{q+1}}{q+1} , \  x\in \mathbb{R}^n.
	\end{align*}
\end{lemma}
\emph{Do similar estimates exist for general HLS type systems \eqref{sys} and \eqref{pde}?} This is also an interesting question, yet as far as we know there has not been any answer to it.

\noindent Proof of comparison principle Lemma \ref{comparisonPrinciple}.

Let $l=(\frac{p+1}{q+1})^{\frac{1}{p+1}}$, $\sigma=\frac{q+1}{p+1}$. So $l^{p+1}\sigma=1$, and $\sigma\leq 1$. Denote $\omega = v- lu^{\sigma}.$
Since
\begin{align*}
	\Delta\omega &= \Delta v - l \nabla\cdot(\sigma u^{\sigma-1}\nabla u)  \\
					&= \Delta v - l\sigma(\sigma-1)|\nabla u|^2 -l\sigma u^{\sigma-1}\Delta u\\
					&\geq -u^q+l\sigma u^{\sigma-1}v^p \\
					&= u^{\sigma-1}((\frac{v}{l})^p - u^{q+1-\sigma}) \\
					&= u^{\sigma-1}((\frac{v}{l})^p - u^{\sigma p}),
\end{align*}
we see $\Delta\omega \geq Cw^p >0$ if $w>0$.
Now, suppose $w>0$ for some $x\in\mathbb{R}^n$.

Case 1: $\exists x_0\in \mathbb{R}^n$ such that $\omega(x_0)=\displaystyle \sup_{\mathbb{R}^n} \omega(x)>0$. Then we have
$\Delta\omega(x_0) \leq 0.$
This contradicts with the fact that $\Delta\omega(x_0) > 0$ since $\omega(x_0) > 0.$

Case 2:There exists a sequence $\{x_R\}$, such that $\displaystyle \lim_{R\rightarrow\infty} \omega(x_R) = \displaystyle \sup_{\mathbb{R}^n} \omega(x) >0$.

Define $\omega_R(x)=\phi(\frac{x}{R})\omega(x)$ where $\phi(x)$ is a cutoff function on $\mathbb{R}^n$ such that $\phi(x)\equiv 1$ on the unit ball $B_1(0)$ and $\phi(x)\equiv 0$ outside $B_2(0)$.  Take $x_R$ such that $\omega_R(x_R) =  \max_{\mathbb{R}^n} \omega_R(x)$.

We see that $\lim_{R\rightarrow\infty} \omega_R(x_R)=\lim_{R\rightarrow\infty} \omega(x_R) = \displaystyle \sup_{\mathbb{R}^n} \omega(x) >0$ and:
\begin{align*}
	& 0=\nabla \omega_R(x_R) = \phi(\frac{x_R}{R})\nabla\omega(x_R) + \frac{1}{R}\nabla\phi(\frac{x_R}{R})\omega(x_R) \\
	\Rightarrow & \nabla\omega(x_R) = O(\frac{1}{R}), \ \text{as } R\rightarrow +\infty.
\end{align*}
So,
\begin{align*}
	0 \geq \Delta\omega_R(x_R)&=\frac{1}{R^2}\Delta\phi(\frac{x_R}{R})\omega(x_R)+\frac{2}{R}\nabla\phi(\frac{x_R}{R})\cdot\nabla\omega(x_R) \\
	& +\phi(\frac{x_R}{R})\Delta\omega(x_R) \\
	\Rightarrow  \Delta\omega(x_R) \leq O(\frac{1}{R^2})
\end{align*}
This contradicts with the relation $\Delta\omega \geq Cw^p >0$ if $w>0$ (suppose $\lim_{R\rightarrow\infty}w(x_R)>C_0>0$ then $\lim_{R\rightarrow\infty}\Delta\omega > C_1>0$ for some $C_1$).

So, we see that $\omega\leq 0$ on $\mathbb{R}^n$, and this finishes the proof. $\Box$

\quad

Lemma \ref{energyEstimates} is first obtained by Serrin and Zou \cite{SZ96} (1996). In \cite{CHL14} a simpler proof is given, to which we refer readers for detail. Here we only sketch the proof.
First, we multiply \eqref{LaneEmden} with $\phi$, the first eigenfunction of $-\Delta$ operator on $B_R$ with eigenvalue $\lambda\sim O(\frac{1}{R^2})$. Then integrate the equations by parts and use the fact that $\phi$ is positive in $B_R$ (which leads to $\frac{\partial \phi}{\partial n}<0$ on $\partial B_R$ by Hopf's Lemma), and we obtain the following inequalities,
\begin{align}\label{uqBoundedByV}
	\int_{B_R} \phi u^q \leq \lambda \int_{B_R} \phi v, \quad
	\int_{B_R} \phi v^p \leq \lambda \int_{B_R} \phi u.
\end{align}
Apply comparison principle lemma \ref{comparisonPrinciple} to the first inequality above to get
	\begin{align*}
		\frac{1}{R^2} \int_{B_R} \phi v\geq C \int_{B_R} \phi v^{\frac{q(p+1)}{q+1}}.
	\end{align*}
Notice that $\frac{q(p+1)}{q+1}>1$ as $pq>1$, so by H\"{o}lder inequality
\begin{align*}
	&\frac{1}{R^2} \int_{B_R} \phi v\geq C (\int_{B_R} \phi v)^{\frac{q(p+1)}{q+1}} R^{-n\frac{qp-1}{q+1}} \\
	\Rightarrow & \int_{B_R} \phi v \leq C R^{n-\beta}.
\end{align*}
Therefore, by \eqref{uqBoundedByV}
\begin{align*}
	\int_{B_R} \phi u^q \leq C R^{n-\beta-2} =CR^{n-q\alpha},
\end{align*}
and \textbf{Case 1: $q \geq 1$}, by H\"{o}lder's inequality:
\begin{align*}
	\int_{B_R} \phi u \leq (\int_{B_R} \phi u^q)^{\frac{1}{q}}(\int_{B_R} \phi )^{\frac{1}{q'}}
			\leq CR^{\frac{n}{q}-\alpha} R^{\frac{n}{q'}}
			= CR^{n-\alpha},
\end{align*}
and by the second inequality of \eqref{uqBoundedByV},
\begin{align*}
	\int_{B_R} \phi v^p \leq CR^{n-\alpha-2} = CR^{n-p\beta}.
\end{align*}
\textbf{Case 2: $q<1$}, this case is more complex than the first one. By the first equation of \eqref{LaneEmden} we have $-\Delta u\leq 0$, then multiply it with $\eta^2 u^{\gamma}$ where $\eta \in C^{\infty}_0(\mathbb{R}^n)$ and $\eta\in(0,1)$ and integrate over whole space, we get
\begin{align*}
	\int_{B_R} \frac{4}{\gamma^2}|D(u^{\frac{\gamma}{2}})|^2\leq \frac{C_\gamma}{R^2}\int_{B_{2R}}u^{\gamma}.
\end{align*}
Now we use Poincar\'{e} inequality to induce an embedding inequality on the support of $\eta$.
Then we use H\"older inequality and estimate for $\int_{B_R} u^q$ in previous proof to obtain,
\begin{align*}
	\int_{B_R} u^{\frac{n}{n-2}\gamma} \leq  C R^{n-\gamma\frac{n}{n-2}\alpha}.
\end{align*}
Then the proof is finished by taking $\gamma\leq q$ then use H\"older inequality to get $$\int_{B_R} u^\theta \leq CR^{n-\theta\alpha,}$$ for any $\theta\in (0,\frac{n}{n-2})$.

\section{Existence and classification of solutions for critical HLS type systems}
The existence of solution to critical HLS type systems, for both integral equations \eqref{sys} and PDE \eqref{pde}, is completely revolved.

In \cite{Lieb83}  E. Lieb (1983) established the existence of ground state solution (i.e. the optimizer of variational problem \eqref{optimizingProblem}). Later, people find that shooting method is a powerful tool to prove existence of solution to both critical and supercritical cases of \eqref{pde} with integer power of Laplacian. We sketch the proof of the shooting method in Section 4.1.

In the following system, Liu, Guo and Zhang \cite{LGZ06} (2006) have obtained the existence for both critical and supercritical cases,
\begin{theorem}\label{LGZ}
The system
\begin{equation} \left \{
 \begin{array}{l}
      (-\Delta)^{k} u = v^q , \;\; u>0, \mbox{ in } R^n,\\
      (-\Delta)^{k} v = u^p, \;\; v>0, \mbox{ in } R^n,
\end{array}
\right.
\end{equation} \label{HLSpde}
admits a positive solution in the critical and super-critical cases
$\frac{1}{p+1}+\frac{1}{q+1}\leq\frac{n-2k}{n}$
for any integer $k$ with $2k<n$. In fact, given $pq>1$, the above system admits a positive radial solution
if and only if $\frac{1}{p+1}+\frac{1}{q+1}\leq\frac{n-2k}{n}$.
\end{theorem}
In Section 4.1, we outline the proof of a more general theorem which contains the result above.
In the scalar case:
\begin{equation}
(-\Delta)^{k} u = u^p,\;\; u>0, \mbox{ in } R^n,
\label{pde10s}
\end{equation}
where $p>1$, $2k<n$. Lei and Li (2013) \cite{LL13} showed that
\begin{theorem}
For $p>1$, equation \eqref{pde10s} admits a locally bounded solution if and only if $p \geq \frac{n+2k}{n-2k}$.
\end{theorem}

\begin{remark}
	Notice that $p>1$ is necessary for the ``only if'' part to be true in the theorem above. If $p=1$ and $k$ is even, \eqref{pde10s} has many solutions. Does \eqref{pde10s} have solution when $k>1$ and $k$ is odd? This is an open question.
\end{remark}
Moreover,
Lei and Li (2013) \cite{LL13} showed that
\begin{theorem} \label{th1.2}
Assume $pq>1$, then the HLS type system (\ref{sys})
has a pair of positive solutions $(u,v) \in L^{p+1}(R^n)
\times L^{q+1}(R^n)$ if and only if it is critical:
\begin{equation} \label{1.3ksce}
\frac{1}{p+1}+
\frac{1}{q+1}=\frac{n-\gamma}{n}.
\end{equation}

\end{theorem}

Now we shall focus on the classification of the solutions. As mentioned in the introduction, for critical HLS type systems, we conjecture that\textbf{ all positive solutions of (\ref{sys}) or (\ref{pde}) are given by a single
family of solutions via translation and scaling.}
As a special case of this conjecture, E. Lieb raised an open problem in \cite{Lieb83} (1983) that asks for $p=q$ if \eqref{sys10} has unique solution up to scaling and translating. Lieb's problem was completely resolved with the introduction of the integral form of method of moving planes by Chen, Li and Ou \cite{CLO06} (2006). The conjecture is still open in general.

The method of moving planes was introduced by A.D. Alexandrov in 1950s, and then developed by J. Serrin \cite{Serrin71} (1971) and Gidas, Ni, Nirenberg \cite{GNN81} (1981). Caffarelli, Gidas and Spruck classified all the solutions to \eqref{pde2} in \cite{CGS89} (1989) by the method of moving planes. Then Chen and Li simplified their proof \cite{CL91} (1991). Wei and Xu generalized this result to higher order conformally invariant equations \cite{WX99} (1999). Then Chen, Li and Ou finally implemented the method of moving planes to the integral equation \eqref{sys} with $p=q$ and solved Lieb's problem \cite{CLO06} (2006).

Hereinafter, by ground state solutions we mean the solution $(u,v)$ corresponding to the optimizer $(f,g)$ of the functional \eqref{optimizingProblem}. To tackle the problem, we break it down to the following open problems:
\begin{enumerate}
\item \emph{all ground state solutions are translation and scaling related;
\item  all finite energy solutions are ground state;
\item  all radial solutions are finite energy solution (the converse is known to be true, see \cite{CLO06});
\item  all solutions decaying to zero at infinity are radial;
\item   all bounded solutions are radial.}
\end{enumerate}

\section{The super-critical HLS type systems}
In this section, we study the existence and asymptotic analysis of solution to supercritical HLS type systems.

\subsection{Existence of solution to critical and supercritical HLS type systems}

For the Lane-Emden system \eqref{LaneEmden}, Serrin and Zou (1998) \cite{SZ98} used shooting method to obtain the existence of solution.
Liu, Guo and Zhang (2006) \cite{LGZ06} introduced a degree approach to shooting method (see also Li (2011) \cite{LiarXiv13}) to obtain radial positive (locally bounded) solution for both critical and supercritical cases of HLS type system \eqref{pde} in a uniform way.

 In short, the degree approach of shooting method combines three ingredients together: degree theory, target map (i.e. shooting) and non-existence on balls with Dirichlet boundary condition (Pohoza\'ev identities). Approach of this kind, relating the existence of solutions in $\mathbb{R}^n$ to the non-existence to a corresponding Dirichlet problem on balls, is implemented by many mathematicians, for instance earlier by Berestycki, Lions and Peletier \cite{BLP81} (1981).

Here we present a theorem by Cheng and Li (2015) \cite{CL15}, which implements degree approach of shooting method to Sch\"odinger type systems with sign-changing nonlinearities (HLS type system \eqref{pde}  is included as a special case). Consider a very general system in whole space,
\begin{align}\label{original}
  \left\{\begin{array}{cl}
      	-\Delta u_{i} = f_{i}(u) & \text{in } \mathbb{R}^n, \\
	       u_{i} > 0             & \text{in } \mathbb{R}^n, \\

        \end{array}
 \right.
\end{align}
and its corresponding local Dirichlet problem,
\begin{align}\label{Dirichlet}
  \left\{\begin{array}{cl}
      	-\Delta u_{i} = f_{i}(u) & \text{in } B_R, \\
	       u_{i} > 0             & \text{in } B_R, \\
           u_i=0                 & \text{on } \partial B_R,

        \end{array}
 \right.
\end{align}
where $B_R=B_R(0)$ for any $R>0$ and $i = 1,2,\cdots, L$. We have
\begin{theorem}\label{existencetheorem}
Given the nonexistence of radial solution to system \eqref{Dirichlet} for all $R>0$, the system \eqref{original} admits a radially symmetric solution of class $C^{2,\alpha}(\mathbb{R}^{n})$ with $0<\alpha<1$, if $f= (f_1(u), f_2(u),\cdots, f_L(u)):\mathbb{R}^L\rightarrow \mathbb{R}^L$ satisfies the following assumptions:
\begin{enumerate}

  \item $f$ is continuous in $\overline{\mathbb{R}^{L}_{+}}$
and locally Lipschitz continuous in $\mathbb{R}^{L}_{+}$, and furthermore,
  \begin{align}\label{decayAssumption}
        \sum_{i=1}^L f_i(u) \geq 0  \text{ in } \mathbb{R}_+^L;
  \end{align}
  \item If $\alpha \in \partial\mathbb{R}^{L}_{+}$ and $\alpha\neq 0$, i.e., for some permutation $(i_1,\cdots,i_L)$, $\alpha_{i_1}=\cdots=\alpha_{i_m}=0$, $\alpha_{i_{m+1}},\cdots, \alpha_{i_L}>0$ where $m$ is an integer in $(0,L)$, then $\exists \delta_0 = \delta_0(\alpha) > 0$ such that for $\beta \in \mathbb{R}^{L}_{+}$ and $|\beta-\alpha|<\delta_0$,
  \begin{align}\label{ControlInequality}
    \sum_{j=m+1}^L |f_{i_j}(\beta)| \leq C(\alpha) \sum_{j=1}^m f_{i_j}(\beta),
  \end{align}
  where $C$ is a non-negative constant that depends only on $\alpha$.
\end{enumerate}
\end{theorem}
Notice that the nonlinear terms of \eqref{HLSpde} (i.e. the HLS type system \eqref{pde} with $\gamma=2k$) satisfy the assumptions \eqref{decayAssumption}-\eqref{ControlInequality}. So Theorem \ref{existencetheorem} recovers Theorem \ref{LGZ} given the non-existence of solution to the corresponding Dirichlet problems. For nonexistence part, people usually apply Pohoza\'ev inequalities \cite{Pohozaev65} (see also \cite{PS86}). This is now more or less standard, a large amount of literature can be found in related topics, for example \cite{Mitidieri93} and \cite{QS07}.

In what follows we outline the proof of Theorem \ref{existencetheorem}, which should exhibit the standard procedure of the degree approach of shooting method.

In the view of seeking radial solutions of \eqref{original}
we solve the following initial value
problem with any initial value $\alpha=(\alpha_1, \cdots, \alpha_L)$ with $\alpha_i>0, i=1,2,....L.$ Denote the solution as $u(r,\alpha)$:
\begin{equation} \label{eq:ode}
\left\{ \begin{aligned}
        &u_i^{''}(r) +\frac{n-1}{r} u_i^{'}(r)= -f_i(u) \\
         &u_i'(0)=0, u_i(0)=\alpha_i \quad i=1,2,\ldots,L.
                          \end{aligned} \right.
\end{equation}
We want to find the suitable initial value $\alpha$ so that $u_i(r,\alpha) >0$ for all $r>0$.

When $L=1$ the question is simple, the assumption that \eqref{Dirichlet} admits no solution
is equivalent to $u_1(r,\alpha) >0$ for all $r>0$. Then there exists a global solution for any initial value, and we are done.

When $L \geq 2$, instead of one dimensional initial value which scales to each other,
we are encountered with multi-dimensional initial value. Among $\alpha_i$'s,
in many critical cases as well as in many supercritical cases,
there is at most one scaling class (one-dimensional) of initial values
from which we can shoot to a global solution. To show the existence of positive solutions of
(\ref{eq:ode}), up to a simple scaling, we have to find the special 1-D initial values.
This is the main reason why there are so many results  in the scalar case
but very little about \eqref{original} for a long time period.

The degree theory approach for the shooting method gives a simple solution to this difficult problem.
It can be used to solve a much larger class of problems.
The argument starts with
defining the target map $\psi$. By $\alpha>0$ we mean $\alpha$ is an interior point of $\mathbb{R}^L_+$, and let $r_0$ be the smallest value of $r$ for which $u_i(r, \alpha)=0$ for some $i$.
Define a map
\begin{align*}
	\psi: & R^L_+ \rightarrow \partial R^L_+	 \\
		& \alpha \mapsto \psi(\alpha)=\left\lbrace \begin{array}{ll}
		u(r_0, \alpha), &\text{ if } \alpha > 0, r_0<\infty,\\
		\displaystyle\lim_{r \rightarrow \infty} u(r, \alpha), &\text{ if } \alpha > 0, r_0=\infty, \\
		\alpha, &\text{ if } \alpha\in \partial R_+^L.
		\end{array}
		\right.
\end{align*}

Then we need to show that $\psi$ is continuous from $R_+^L$
to $ \partial R_+^L$.
Assumptions \eqref{decayAssumption}-\eqref{ControlInequality} guarantee this.

In the next step, applying the degree
theory, we show that $\psi$ is onto from $A_a$ to $B_a$ where:

\begin{equation} \label{eq:scalar1}
\left\{ \begin{aligned}
 & A_a\triangleq \{ \alpha \in R_+^L \,\, {\bf \mid } \quad {\tiny \displaystyle\sum_{i=1,\cdots, L}} \alpha_i=a\},\\
  & B_a\triangleq\{ \alpha \in \partial R_+^L \,\, {\bf \mid } \quad {\tiny \displaystyle\sum_{i=1,\cdots, L}} \alpha_i \leq a\},
                          \end{aligned} \right.
\end{equation}
for any $a>0$.
In particular, there exists at least one
$\alpha_a \in A_a$ for
sufficiently small $a>0$ such that $\psi(\alpha_a)=0$.

Shooting from the initial value $\alpha_a$, by the assumption that the system
\eqref{Dirichlet} admits no radially symmetric solution, we obtain a solution of \eqref{original} which decays to 0 at infinity ($\displaystyle\lim_{r \rightarrow \infty} u(r, \alpha_a)=\psi(\alpha_a)=0$). Notice that we did not rule out the radial solution that does not decay at infinity for system with sign-changing nonlinearity, however for positive nonlinear source term such as HLS type system \eqref{pde}, the radial solution must monotone decrease, hence its radial solution must decay to zero at infinity.

We remark that assumption \eqref{ControlInequality} is a necessary condition for the target map to be continuous. A trivial degenerate case in the following does not satisfy \eqref{ControlInequality},
\[ \left\{ \begin{array}{l}
         -\triangle u_1 = u_1^p, \\
                  -\triangle u_2=u_2^q,
                  \end{array} \right.
         \]
which can be decoupled to the study of two scalar equations.
If $p> \frac{n+2}{n-2}$ and $q< \frac{n+2}{n-2}$ then the target map is not continuous at
$(0,a)$ for any $a>0$. In fact, $\phi(0,a)=(0,a)$ and $\psi(\delta, a)=(h(\delta),0)$ where $0<h(\delta)<\delta$
and $\psi(\delta, a)\rightarrow (0,0) \neq (0,a)$ as $ \delta \rightarrow 0^+$.

\subsection{Asymptotic analysis}

In the asymptotic analysis we try to answer two questions: What are the asymptotic expansions of the radial solutions? Are all radial solutions scaling related?

In \cite{LLM12}, Lei, Li and Ma (2012) considered a more general system as following,
\begin{equation}
 \left \{
   \begin{array}{l}
      u(x) = \frac{1}{|x|^{\alpha}}\int_{R^{n}} \frac{v(y)^q}{|y|^{\beta}|x-y|^{\lambda}}  dy\\
      v(x) = \frac{1}{|x|^{\beta}}\int_{R^{n}} \frac{u(y)^p}{|y|^{\alpha}|x-y|^{\lambda}}  dy
   \end{array}
   \right.\label{ELWHLS}
 \end{equation} where $\alpha+\beta+\lambda \leq n$, and
     \begin{equation}
     \left \{
   \begin{array}{l}
   u, v \geq 0, \;\; 0 <p, q<\infty,\;
    0 <\lambda < n, \; \alpha+\beta  \geq 0,\; \\
   \frac{\alpha}{n}  <
  \frac{1}{p+1}<\frac{\lambda+\alpha}{n},\;\frac{1}{p+1}+\frac{1}{q+1}=\frac{\lambda+\alpha+\beta}{n}.
\end{array}
   \right.
 \label{ELWHLScond3}
   \end{equation}
Notice that HLS type system \eqref{sys} is a special case ($\alpha=\beta=0$) of \eqref{ELWHLS}.

\begin{definition}
 A function $u$ is asymptotic to $\frac {B}{|x|^t}$ near
$x=\infty$, that is
\begin{equation}
u(x) \simeq \frac {B}{|x|^t}\;{\rm at} \;|x| \simeq \infty, \quad
{\rm if} \lim_{|x|\to \infty}|x|^t u(x) = B,\label{def_2}
\end{equation} for a positive number $t$ and $B \neq 0\; (\infty)$.
\end{definition}

Lei, Li and Ma proved that

\begin{theorem}
\it Let $(u,v) \in L^{p+1}(R^n)\times
L^{q+1}(R^n)$ be a pair of positive solutions of the system
(\ref{ELWHLS}) with (\ref{ELWHLScond3}). Suppose that $p \geq 1,\;
q \geq 1,\; pq \neq 1$ and $\alpha+\beta \geq 0$. If $\lambda
q+\beta (q+1)>n$, then  we have
\begin{equation}
     u(x) \simeq  \frac {B_0}{|x|^{\lambda+\alpha}}, \text{ for large $|x|$,}\\
     \label{uasyminf}
\end{equation}

and

\begin{eqnarray}
     v(x) & \simeq & \left\{\begin{array}{cll}
          \displaystyle\frac {B_1}{|x|^{\lambda+\beta}},
              \quad   &{ if}\; \lambda p +\alpha(p+1)>n\\
        \displaystyle\frac {B_2 |\ln |x||}{|x|^{\lambda+\beta}},
             \quad &{ if}\; \lambda p +\alpha(p+1)=n \\
          \displaystyle\frac {B_3} {|x|^{(\alpha +\lambda)(p+1)+\beta-n}},
              \quad &  { if}\; \lambda p +\alpha(p+1)<n   \label{vasyminf}
    \end{array}\right. \text{ for large $|x|$}
 \end{eqnarray}
where $B_0 =\int_{R^n}\frac{v^q(y)} {|y|^{\beta}} \,dy$, $B_1
=\int_{R^n}\frac{u^p(y)} {|y|^{\alpha}} \,dy$, $B_2 =
|S^{n-1}|(\int_{R^n}\frac{v^q(y)} {|y|^{\beta}} \,dy)^p$ and\\
$B_3 = (\int_{R^n}\frac{v^q(y)} {|y|^{\beta}} \,dy)^p
\int_{R^n}\frac{dz}{|z|^{2n-(\alpha+\lambda)(p+1)}|e-z|^{\lambda}}$.

\end{theorem}
\par Recently, the authors consider the supercritical Lane-Emden system \eqref{LaneEmden}
with $p,q$ satisfying the supercritical condition, i.e.
\begin{equation}\label{4.2.2}
\frac 1{p+1}+\frac 1{q+1}<1-\frac2{n}, \quad 0<p,q<+\infty.
\end{equation}
The existence of the positive radial symmetric solutions are already known by \cite{LGZ06} via shooting method. We consider the asymptotic property of these radial solutions in the supercritical cases.
\begin{theorem}
\label{thm4.2.1}
Suppose $(u,v)$ be locally bounded non-negative radial solutions of \eqref{LaneEmden} with $p,q$ satisfying \eqref{4.2.2}. Then either
\begin{equation}\label{4.2.13}
C_1 r^{-\frac{2(p+1)}{pq-1}}\leq u(r)\leq C_2 r^{-\frac{2(p+1)}{pq-1}},\quad C_1 r^{-\frac{2(q+1)}{pq-1}}\leq v(r) \leq C_2 r^{-\frac{2(q+1)}{pq-1}}
\end{equation}
for some positive constants $C_1,C_2>0$ and $r$ large or $u=v\equiv 0$.
\end{theorem}
A more interesting open problem is if $u(r) r^{\frac{2(p+1)}{pq-1}}, v(r)r^{\frac{2(q+1)}{pq-1}}$ converge as $r\rightarrow\infty$.

Proof.
Recall that we define $\alpha=\frac{2(p+1)}{pq-1},\beta=\frac{2(q+1)}{pq-1}$ (see Section 2).
Considering the radial solutions of \eqref{LaneEmden}, we have the following ODE system,
\begin{align}\label{odeLaneEmden}
\left\lbrace \begin{array}{cc}
 -(r^{n-1}u')'=r^{n-1}v^p, \\
 -(r^{n-1}v')'=r^{n-1}u^q,
\end{array}\right.
\end{align}
Since we consider the radial solutions, by Lemma \ref{energyEstimates}, the upper bound of \eqref{4.2.13} follows easily.

\par It remains to prove that $u,v$ are positive solutions, then $\displaystyle\liminf_{R\rightarrow +\infty} R^\alpha u(R)=c>0$.
\par Suppose not, we may assume $u(r)=w(r)r^{-\alpha}$ and $\displaystyle\liminf_{r\rightarrow +\infty} w(r)=0$. Now, we derive the upper bound of $-u'(r), -v'(r)$. Again, by comparison principle, we have
$$v(r)\leq \sigma u^{\frac{q+1}{p+1}}=\sigma w^{\frac{q+1}{p+1}}(r)r^{-\beta},$$
where $\sigma=\sqrt[p+1]{\frac{p+1}{q+1}}$. As discussed in \cite{Mitidieri92}, we have
\begin{equation}
\begin{split}
-ru'(r)=(N-2)u(r)-(N-2)u(0)-\int_r^\infty sv^p(s)ds\leq (N-2)u(r)=(N-2)w(r)r^{-\alpha}
\end{split}
\end{equation}
The same conclusion also holds for $-rv'(r)$ such that $-rv'(r)\leq (N-2)\sigma w^{\frac{q+1}{p+1}}(r) r^{-\beta}$.
Recall the following Pohozaev's identity for radial solutions,
\begin{equation}
\begin{split}
&(\frac n{p+1}-a_1)\int_0^R s^{n-1} v^{p+1}(s)ds+(\frac{n}{q+1}-a_2)\int_0^R s^{n-1}u^{q+1}(s)ds\\=
&R^{n-1}(\frac{Rv^{p+1}(R)}{p+1}+\frac{Ru^{q+1}(R)}{q+1}+Ru'v'(R)+a_1u'v(R)+a_2v'u(R)).
\end{split}
\end{equation}
We distinguish two cases:
\par 1. There exists $r_0>0$ such that for $r>r_0$, $w'(r)\leq 0$.
\par 2. For any $l>0$, there exists $s_l>l$ such that $w'(s_l)>0$. Since $\displaystyle\liminf_{r\rightarrow \infty} w(r)=0$,
we can get a sequence of local minimum points $r_l$ of $w$, such that $w(r_l)\rightarrow 0$ as $r_l\rightarrow +\infty.$
\par For both cases, we have
\begin{equation}\label{lowerboundAsymptotic}
-u'(r_l)=-w'(r_l)r^{-\alpha}_l+\alpha w(r_l)r_l^{-\alpha-1}\geq \alpha w(r_l)r_l^{-\alpha-1}.
\end{equation}
Now consider Pohozaev's identity at $r_l$ and $a_1=\frac n{p+1}$. By the supercritical property of $p,q$, we have
\begin{equation}
\begin{split}
\int_0^{r_l} r^{n-1}u^{q+1}(r)dr \leq C(w^{q+1}(r_l)+w^{1+\frac{q+1}{p+1}}(r_l))r_l^{n-2-\alpha-\beta}\leq Cw^{1+\frac{q+1}{p+1}}(r_l)r_l^{n-2-\alpha-\beta}
\end{split}
\end{equation}
where we have used $\frac{q+1}{p+1}<q$.  By comparison principle, we also have
$$\int_0^{r_l} r^{n-1}v^{p+1}(r)dr \leq Cw^{1+\frac{q+1}{p+1}}(r_l)r_l^{n-2-\alpha-\beta}.$$
 Now we return to the first equation of \eqref{odeLaneEmden},
\begin{equation}\label{upperboundAsymptotic1}
\begin{split}
-r_l^{n-1}u'(r_l)&=\int_0^{r_l}s^{n-1} v^p(s)ds\\
&\leq \left(\int_0^{r_l} s^{n-1}v^{p+1}\right)^{\frac p{p+1}}\left(\int_0^{r_l} s^{n-1}ds\right)^{\frac 1{p+1}}\\
&\leq C w^{\frac p{p+1}+\frac{p(q+1)}{(p+1)^2}}(r_l)r_l^{n-2-\alpha}
\end{split}
\end{equation}
Combine \eqref{lowerboundAsymptotic} and \eqref{upperboundAsymptotic1} we get $$C w^{\frac p{p+1}+\frac{p(q+1)}{(p+1)^2}}(r_l)r_l^{n-2-\alpha} \geq -r_l^{n-1}u'(r_l)\geq \alpha w(r_l) r_l^{n-2-\alpha}.$$
From $w(r_l)\neq 0$, we obtain
\begin{equation}
w^{\frac{pq-1}{(p+1)^2}}(r_l)\geq \frac 1C
\end{equation}
This contradicts with $w(r_l)\rightarrow 0$.
Hence $u(r)\sim r^{-\alpha}$ near $\infty$.
Then integrate the second equation of \eqref{odeLaneEmden} from 0 to $r$,
\begin{equation}
-r^{n-1} v'(r)=\int_0^r s^{n-1} u^q ds\geq Cr^{n-q\alpha}.
\end{equation}
Dividing the above inequality by $r^{n-1}$ and integrating from $r$ to $\infty$, we derive $v(r)\geq c r^{-\beta}$.
\par This ends the proof of the asymptotics of $u,v$ at $\infty$.
$\Box$

\section{Subcritical HLS type systems: Liouville type theorem and the Lane-Emden conjecture}

For subcritical HLS type systems, various Liouville type theorems are obtained. This kind of results are often based on the study of the Lane-Emden system \eqref{LaneEmden}.
The Lane-Emden conjecture has been lasting unsolved for decades. Many mathematicians have contributed in this question, for example, the pioneer job done by Mitidieri (1992) \cite{Mitidieri92} (see also \cite{Mitidieri96}) which solves the Lane-Emden conjecture in radial case. For expository reference about the Lane-Emden conjecture, readers can check \cite{CHL14, Souplet09} and reference therein. Among these mathematical works we mention a couple of results below.

	For $pq\leq 1$, it is known that system \eqref{LaneEmden} has no positive classical supersolutions (see Serrin and Zou (1996) \cite{SZ96}).

For $n=3$,  the conjecture is solved by two papers. First, Serrin and Zou (1996) \cite{SZ96} proved that there is no positive solution with polynomial growth at infinity.
\begin{theorem}[Serrin-Zou-1996]\label{SZ06}
		Let $n=3$. Lane-Emden system \eqref{LaneEmden} admits no solution given the solution has at most polynomial growth at infinity.
	\end{theorem}
Then Pol\'{a}\v{c}ik, Quittner and Souplet (2007) \cite{PQS07} removed the growth condition. In fact, they proved that no bounded positive solution implies no positive solution.
\begin{theorem}[Pol\'{a}\v{c}ik-Quittner-Souplet-2007 ]\label{PQS09}
		Let $pq>1$. Assume that \eqref{LaneEmden} does not admit
		any bounded nontrivial (nonnegative) solution in $\mathbb{R}^n$.
		Let $\Omega\neq\mathbb{R}^n$ be a domain. There exists $C=C(n,p,q)>0$
		such that any solution $(u,v)$ of \eqref{LaneEmden} in $\Omega$
		satisfies
		\begin{align*}
			u(x)\leq C dist^{-\alpha}(x,\partial\Omega), \ x\in\Omega,
		\end{align*}
		and
		\begin{align*}
			v(x)\leq C dist^{-\beta}(x,\partial\Omega), \ x\in\Omega.
		\end{align*}
	\end{theorem}

	\begin{remark}
	In \cite{PQS07}, $p,q$ were assumed to be both $>1$, however, their proof is valid for $pq>1$ and can be verified directly.
	\end{remark}

This result has two important consequences. One is that combining with Serrin and Zou's result, one can prove the conjecture for $n=3$.
	\begin{cor}
	Lane-Emden conjecture is true when $n=3$.
	\end{cor}
The other consequence is that one can reduce the Lane-Emden conjecture to the nonexistence of bounded positive solution.
Thus, hereinafter we always assume that $(u,v)$ are bounded.

For $n=4$, the conjecture is recently solved by Souplet (2009) \cite{Souplet09}. In \cite{SZ96}, Serrin and Zou used the integral estimates to derive the nonexistence results. Souplet further developed the approach of integral estimates and solved the conjecture for $n=4$ along the case $n=3$. In higher dimensions, this approach provides a new subregion where the conjecture holds, but the problem of full range in high dimensional space still seems stubborn.
\begin{theorem}[Souplet-2009]
	For subcritical Lane-Emden system \eqref{LaneEmden}, i.e. $p,q>0$ and $\frac{1}{p+1} +\frac{1}{q+1} > \frac{n-2}{n}$,
	\begin{enumerate}
	\item  $n= 3$ or 4, the system \eqref{LaneEmden} has no positive classical solutions.

	\item  $n\geq 5$, $pq>1$,  along with
	\begin{align}
		\max\{\alpha,\beta\} > n-3,
		\end{align}
	then system \eqref{LaneEmden} has no positive classical solutions.
	\end{enumerate}

	\end{theorem}
For higher order HLS type systems \eqref{pde}, Arthur, Yan and Zhao (2014) \cite{AYZ14} have proved a Liouville theorem with an adapted idea of measure
and feedback argument developed by Souplet (2009) \cite{Souplet09}.

Existence of solutions to subcritical Lane-Emden system with double bounded coefficients can be extended from low dimension to system in higher dimension which fails the integral estimates \eqref{energyEstimates}. This may imply that, integral estimates are essential to problem of nonexistence of solutions. In \cite{CHL14}, Cheng, Huang and Li raised the conjecture below, which is proven to be equivalent to the Lane-Emden conjecture.
\begin{conjecture}
For solution $(u,v)$ to the Lane-Emden system with $p\geq q$, there exist an $s>0$ such that $n-s\beta<1$ and
\begin{align}
	\int_{B_R} v^s \leq C R^{n-s\beta}.  \label{vEstimateBetter}
\end{align}
\end{conjecture}

	\begin{theorem}[Cheng-Huang-Li-2014]\label{CHL14theorem}
		Let $(u,v)$ be a non-negative bounded solution to \eqref{LaneEmden}, with
		$\frac{1}{p+1} +\frac{1}{q+1} > \frac{n-2}{n}$.
		Assume additionally $p\geq q$ and there exists an $s>0$ satisfying $n-s\beta<1$ such that v satisfies \eqref{vEstimateBetter}
		then $u, v \equiv 0$.
	\end{theorem}

\begin{remark}
Energy estimate \eqref{vEstimateBetter} is a necessary condition to the Lane-Emden conjecture. One just needs to notice that when $u,v\equiv 0$, \eqref{vEstimateBetter} is obviously satisfied. The key to the proof of Theorem \ref{CHL14theorem} is to show \eqref{vEstimateBetter} is sufficient.
\end{remark}

\begin{remark}
The assumption on $v$ is weaker than the corresponding assumption on $u$ due to a comparison principle between $u$ and $v$.

Indeed, we can  replace \eqref{vEstimateBetter} by: for some $r>0$, $n-r\alpha <1$,
\begin{align}
	\int_{B_R} u^r \leq CR^{n-r\alpha}.  \label{uEstimateBetter}
\end{align}
\end{remark}

The proof of Theorem \ref{CHL14theorem} is also based on the feedback argument developed by Souplet, though some basic estimates are adapted as needed. Here a
sketch of the proof of Theorem \ref{CHL14theorem} is presented.
\begin{enumerate}
\item Let $F(R) = \int_{B_R} u^{q+1}.$
By a Rellich-Pohozaev identity of \eqref{LaneEmden} on $B_R$, one can estimate $F(R)$ by quantities on $S^{n-1}$,
\begin{align*}
	F(R) \leq G_1(R) + G_2(R),
\end{align*}
where
\begin{align}
	& G_1(R)=R^n\int_{S^{n-1}}u^{q+1}(R), \label{g1}  \\
	    & G_2(R)=R^n\int_{S^{n-1}}\left( |D_x u(R)|+R^{-1}u(R)\right) \left( |D_x v(R)|+R^{-1}v(R)\right) .  \label{g2}
\end{align}

\item Heuristically, we prove that there exists a sequence $\{R_j\} \rightarrow \infty$ such that the following estimate holds,
\begin{align*}
	F(R_j) \leq C R_j^{-a}F^b(R_j),
\end{align*}
with $a>0$ and $b<1$. Then $F(R)\equiv 0$.
\end{enumerate}

	We start with estimate on $G_1(R)$. By H\"{o}lder's inequality and Sobolev embedding on $S^{n-1}$,
	\begin{align}\label{holderInequalityG1}
	    \|u\|_{q+1} &\leq \|u\|_{\lambda}^{\theta}\|u\|_{\mu}^{1-\theta} \\
	                &\leq (R^2\|D_x^2 u\|_l+\|u\|_1)^{\theta}(R^2\|D_x^2 u\|_k+\|u\|_1)^{1-\theta},
	\end{align}
	So,
	\begin{align*}
		G_1(R) \leq R^n \left( (R^2\|D_x^2 u\|_l+\|u\|_1)^{\theta}(R^2\|D_x^2 u\|_k+\|u\|_1)^{1-\theta}\right) ^{q+1}.
	\end{align*}
	Then by $W^{2,p}$-estimate, energy estimates in Lemma \ref{energyEstimates} together with the assumed integral estimate we have, $\exists \tilde{R}\in(R,2R)$ such that
	\begin{align*}
	    G_1(\tilde{R})\leq R^{-a}F^b(4R).
	\end{align*}

	\begin{remark}
	The existence of such $\tilde{R}$  is guaranteed by a fact that for $f\in L_{loc}^p(B_R)$, $\exists \tilde{R}\in(R,2R)$ such that
	\begin{align*}
		\|f\|_{L^{p}(S^{n-1})}(\tilde{R})  \leq M R^{-\frac{n}{p}}\|f\|_{L^{p}(B_{2R})}.
	\end{align*}
	Notice that the same conclusion can be made to finitely many $L^{p_i}_{loc}$ functions.
	\end{remark}

Moreover, we have the following relation between $a$ and $b$,
	\begin{align*}
		  a = (\alpha+\beta+2-n)(1-b).
	\end{align*}
So, to show that $a<0$ and $b<1$ we only need to verify that $b<1$, and this is guaranteed by
	\begin{align*}
		b<1 \Leftrightarrow n-1 < s\beta.
	\end{align*}
Then we look at estimate on $G_2$.
By H\"{o}lder's inequality, \eqref{g2} becomes,
\begin{align*}
	G_2(R) &\leq R^n\||D_x u| + R^{-1}u\|_z\||D_x v| + R^{-1}v\|_{z'} \\
				&\leq R^n(\|D_x u\|_z+R^{-1}\|u\|_1)(\|D_x v\|_{z'} + R^{-1}\|v\|_1),
\end{align*}
Again applying  H\"{o}lder's inequality and Sobolev embedding to $\|D_x u\|_z$ and $\|D_x v\|_{z'}$, we get
\begin{align*}
    G_2(R)\leq & C R^{n+2}(\|D_x^2 u\|_{l} + R^{-1}\|D_x u\|_1 + R^{-2}\|u\|_1)^{\tau_1} \\
            & \times(\|D_x^2 u\|_k + R^{-1}\|D_x u\|_1 + R^{-2}\|u\|_1)^{1-\tau_1} \\
            & \times(\|D_x^2 v\|_{1+\epsilon} + R^{-1}\|D_x v\|_1 + R^{-2}\|v\|_1)^{\tau_2} \\
            & \times(\|D_x^2 v\|_m + R^{-1}\|D_x v\|_1 + R^{-2}\|v\|_1)^{1-\tau_2},
\end{align*}
with
\begin{align}\label{rangeZ}
	\max\left\lbrace \frac{1}{k} - \frac{1}{n-1}, \frac{1}{n-1}\right\rbrace \leq \frac{1}{z} \leq \min \left\lbrace \frac{1}{l}-\frac{1}{n-1}, \frac{1}{q+1}+\frac{1}{n-1}\right\rbrace.
\end{align}
Similar to $G_1(R)$,
by $W^{2,p}$-estimate, basic integral estimates \eqref{uvEstimate1} and \eqref{uvEsitmatePQ} together with the assumed integral estimate we have, $\exists \tilde{R}\in(R,2R)$ such that
	\begin{align*}
	    G_2(\tilde{R})\leq R^{-\tilde{a}}F^{\tilde{b}}(4R).
	\end{align*}
(Surprisingly) similar to $a,b$ for $G_1(R)$, $\tilde{a}$ and $\tilde{b}$ also have the relation
	\begin{align*}
		\tilde{a} = (\alpha+\beta+2-n)(1-\tilde{b}),
	\end{align*}
	so $\tilde{a}>0$ if and only if $\tilde{b}<1$,
	and
	\begin{align}
	\tilde{b} &< 1 \label{bLessThan1} \\
	\Leftrightarrow & (m-\frac{k}{l})\frac{1}{z}+(\frac{k}{n-1}+(m-1)(k-1))\frac{1}{l}+\frac{m-2}{n-1}-(m-1)>0, \label{zAndl}
	\end{align}
	Now we only need to verify this inequality holds when $\frac{1}{z}$ takes value as its upper bounds in \eqref{rangeZ}:
	\begin{align*}
		\frac{1}{z} \leq \min \left\lbrace \frac{1}{l}-\frac{1}{n-1}, \frac{1}{q+1}+\frac{1}{n-1}\right\rbrace.
	\end{align*}

\textbf{Case 1.} If $\frac{1}{l}-\frac{1}{n-1} \geq \frac{1}{q+1}+\frac{1}{n-1}$,
\begin{align*}
 \tilde{b} < 1
 	\Leftrightarrow  n- s\beta  < 1,
\end{align*}
which is true under our assumption.

\textbf{Case 2.} If $\frac{1}{l}-\frac{1}{n-1} < \frac{1}{q+1}+\frac{1}{n-1}$,
\begin{align*}
 \tilde{b} < 1
 	\Leftrightarrow  \frac{1}{k} < \frac{1}{l} < \frac{2}{n-1}+\frac{1}{q},
\end{align*}
where the inequalities on the right are true since $\frac{1}{l} < \frac{1}{q+1}+\frac{2}{n-1}< \frac{2}{n-1}+\frac{1}{q}$.

 In all, \eqref{bLessThan1} always holds under our assumption $n-s\beta<1$. This finishes the proof of Theorem \ref{CHL14theorem}.

\newpage

\providecommand{\bysame}{\leavevmode\hbox to3em{\hrulefill}\thinspace}
\providecommand{\MR}{\relax\ifhmode\unskip\space\fi MR }
\providecommand{\MRhref}[2]{%
  \href{http://www.ams.org/mathscinet-getitem?mr=#1}{#2}
}
\providecommand{\href}[2]{#2}


\section{Introduction}
There have been many algorithmic approaches proposed to solve complex optimisation problems, including constrained optimisation problems (COP). Several approaches have been proposed to tackle the constraints in constrained problems. Most of the research has been focused on introducing differential evolution (DE) \cite{storn1997differential}, particle swarm optimisation (PSO) \cite{eberhart1995new} and evolutionary strategies (ES) \cite{schwefel1993evolution} to solve numerical optimisation problems. In order to deal with these constrained problems, there have been techniques that applied to these algorithms such as penalty functions, special operators (separating the constraint and objective function treatment) and decoder based methods. We refer the reader for a survey of constraint handling techniques in evolutionary computing methods to \cite{mezura2011constraint}.

In order to compare and evaluate the evolutionary algorithms many approaches have been used. One is finding which algorithm performs better on a set of continuous problems using benchmarks sets \cite{mallipeddi2010problem,hansen2010real}. Recently, there has been an increasing interest to analyse the problem features that make it hard to solve. Initial studies have been carried out in the field of continuous optimisation in \cite{mersmann2010benchmarking}. Furthermore, there have been techniques that generate a variation of problem instances from easy to hard. Then, the features of this problem instances are analysed in order to find which of them make the problems hard or easy to solve. Generating the variety of problem instances from easy to hard ensures that the knowledge obtained from analysis is reliable.

Although there is not only a standalone feature that makes a problem hard to solve, but it is assumed that constraints are very important in constrained continuous problems. The evolving approach that has been used to analyse the constraint features and their effects on COP's difficulty is discussed in \cite{poursoltan2014feature,shayan2015}. The idea is to evolve constrained problem instances (by using an evolutionary algorithm) in order to identify the constraint features with more contribution to problem difficulty.

In this paper, by using a single-objective evolutionary algorithm, we generate hard and easy COP instances for DE, ES and PSO algorithms. Later, we solve the generated instances using one algorithm by the other algorithms. The results show that the hardest generated instances using one algorithm are still hard for the other ones. To get better insight, we use multi-objective evolving approach to generate instances that are hard for one algorithm but still easy for the others. By analysing how an algorithm fails in conditions where the rest perform well, we can derive its strengths and weaknesses over constraint features. Our study shows the effectiveness of constraint features that make the problems hard for one and easy for the other algorithms. It can be translated as over which features of constraints, they make the problems hard for a certain algorithm but still easy for the others.

The remainder of this paper is as follows: In Section 2 we introduce the concept of COPs. Then we discuss the evolver (single and multi-objective evolutionary approach) and the solver algorithms (DE, ES and PSO) we use in our experiments. In Section 3 we analyse the performance of various algorithms on each others hard and easy instances (using the single-objective evolver). Section 4 includes the multi-objective approach that generates hard instances for one but easy for the other algorithms. Furthermore, we carry out the analysis of linear and quadratic constraint features that make the problem hard for one and still easy for the rest. Finally, we conclude with some remarks.

\section{Preliminaries}
\subsection{Constrained continuous optimisation problems}
In this study, constrained continuous optimisation problems with inequality and equality constraints are investigated. These problems are optimisation problems where a function $f(x)$ should be optimised with respect to a given set of constraints.

Single-objective functions $f \colon S \rightarrow \mathbb{R}$ with $S \subseteq \mathbb{R}^n$ are considered in this research.
The constraints impose a feasible subset $F \subseteq S$ of the search space $S$ and the aim is finding $x \in S \cap F$ which minimises $f$. Formally, we state the problems as follows:

	\begin{equation}
		\begin{split}
		\text{minimize } &\quad    f(x), \quad	x = (x_1,\ldots,x_n) \in \mathbb{R}^n  \\
		\text{subject to}&\quad	g_i(x)  \leq 0 \quad  \forall   i \in \{1,\ldots,q\}\\
		&	\quad	h_j(x) = 0 \quad \forall j \in \{q + 1,\ldots, p\}
		\label{eq:f}
		\end{split}
		\end{equation}

where $x = (x_1,x_2,\dots,x_n)$ is an $n$ dimensional vector and $x \in S  \cap F$. The $g_{i}(x)$ (inequality) and $h_{j}(x)$ (equality) constraints could be linear/nonlinear. Also,
the equality constraints are usually replaced by $|h_j(x)|  \leq \epsilon$ where $\epsilon=10e^{-4}$ \cite{mallipeddi2010problem}. The feasible region $F \subseteq S$ of the search space $S$ is defined by

\begin{equation}
l_i  \leq x_i  \leq u_i,\quad  \quad     1 \leq i  \leq n
\label{equality}
\end{equation}

\noindent where $l_i$  and $u_i$ denote lower and upper bounds respectively for the $i$th variable in which  $1\leq i \leq n$.
In this paper, we focus on the ability of constraints (linear, quadratic) to make a problem hard or easy. The features of these constraints and their effect on problem difficulty is discussed. The constraints are of the following form:

\begin{equation}
\text{linear constraint} \quad	g(x) = b + a_{1}x_{1} + \ldots + a_{n}x_n
\label{lineareq}
\end{equation}

\begin{equation}
\text{quadratic constraint} \quad	g(x) = b +a_{1}x_{1}^{2} + a_{2}x_{1}\ldots + a_{2n-1}x_n^{2} +a_{2n}x_n
\label{quadraticeq}
\end{equation}

\noindent or a combination of them, where $x_{1},x_{2} \dots ,x_{n}$ are values from Equation 1 and $a_{1},a_{2},\dots,a_{n}$ are coefficients within lower ($l_{i}$) and upper bounds ($u_{i}$). We assume univariant quadratic function to analyse each $x_{n}$ (with exponent 2) independently. Also, unvivarient quadratic constraints are more popular in recent benchmarks \cite{mallipeddi2010problem}.
In order to include the optimum of objective function in feasible area, we set $b\leq0$ (we assume the objective function optimum is zero).
\subsection{Algorithms}
We now introduce the algorithms for constrained continuous optimisation that are subject to our investigation.

One of the most prominent evolutionary algorithms for COPs is $\epsilon$-constrained differential evolution with an archive and gradient-based mutation ($\epsilon$DEag). The algorithm is the winner of 2010 CEC competition for continuous COPs \cite{mallipeddi2010problem}. The $\epsilon$DEag uses $\epsilon$-constrained method to transform algorithms for unconstrained problems to constrained ones. It adopts $\epsilon$-level comparison to order the possible solutions. In other words, the lexicographic order is used in which constraint violation ($\phi(x)$) has more priority and proceeds the function value ($f(x)$). For more details we refer the reader to \cite{takahama2010constrained}.

The second algorithm we use in this paper is a $(1+1)$ CMA-ES for constrained optimisation \cite{arnold20121+}. The $(1+1)$ CMA-ES in \cite{igel2006computational} is a variant of  $(1+1)$-ES which adapts the covariance matrix of its offspring distribution in addition to its global step size. The idea behind the constraint handling approach of this algorithm is to obtain approximations to the normal vectors directions in the vicinity of the current solutions locations by low-pass filtering steps which violates the respective constraints and reducing the variance of the offspring distribution in these directions. Incorporating this constraint handling approach with $(1+1)$ CMA-ES makes an algorithm which is significantly more efficient than other approaches for constrained evolutionary algorithms. Also, the selected algorithm is not sensitive to the rotation of the problem search space. We refer the reader to \cite{arnold20121+} for more details and implementation.

The third algorithm that is used in our investigation is a particle swarm optimisation. This algorithm (HMPSO) applies a method that uses parallel search operator in which it divides the current swarm into various sub-swarms and locates the solution between them. In each sub-swarm, all particles follow the local best (fittest particle) which improves them to be more fitter. Also, since all sub-swarms are located around different optima (in parallel), then it is more possible to locate multiple optima which improves the diversity of algorithm. Dividing the swarms into sub-swarms improves the diversity of the algorithm. Also, choosing the local best in each sub-swarm can attract the other particles to fitter positions. We refer the reader to \cite{wang2009hybrid} for detailed algorithm and implementation.

\subsection{Features of Constraints}
In this paper we analyse the constraint features of generated problem instances. These features are constraint coefficients relationships such as standard derivation, angle between constraint hyperplanes, feasibility ratio in vicinity of optimum, number of constraints, shortest distance of constraint hyperplane to optimum. The details of these features are discussed in \cite{shayan2015}.

\section{Single-objective Investigations}

We first consider different algorithms and compare their relative performance on each other's generated hard and easy instances. We use single-objective evolver to evolve and generate hard and easy instances for all types of algorithms. The detailed procedure and results for DE instances are discussed in \cite{shayan2015}. For this experiment, we perform 30 independent runs generating easy and hard instances for PSO and ES solvers. It means, the single-objective evolver only generates instances that are hard/easy for one type of algorithm (PSO, ES and DE). The required function evaluation number (FEN) for solving these instances (PSO, ES and DE) is used as fitness value for single-objective evolver. The parameters for solvers are identical to \cite{takahama2010constrained,wang2009hybrid,arnold20121+}. Also, we run our experiments on Sphere function (bowl shaped)\cite{hansen2010real}. We now have three groups of easy and hard instances generated for DE, ES and PSO algorithms. We then compare the DE, ES and PSO algorithms by applying them on each other's easy and hard instances. The analysis is done by comparing the required FEN for an algorithm to solve the other's generated problem instances. Then, it is possible to derive strengths and weaknesses of the considered algorithms by observing how well one algorithm performs in conditions where the other algorithms fail (or it is difficult for them).
Table \ref{table:perflinearsphere} and \ref{table:perfquadsphere} show different algorithms performance on Sphere objective functions with linear/quadratic constraints (1 to 5 constraints). We also run our experiments on different objective functions such as Ackley and Rosenbrock. The results are shown in Tables \ref{table:perflinearackley}, \ref{table:perfquadackley}, \ref{table:perflinearrosenbrock} and \ref{table:perfquadrosenbrock}. It is interesting that all objective functions follow similar pattern.  Considering the required FEN to solve each instances, it is observed that hard instances are still the hardest for their own algorithms and hard for the others. It implies that the hard instances share some common features to make it difficult to solve for all solvers. However, the obtained knowledge is not enough to compare the algorithm capabilities to solve hard problem instances.

  \section{Multi-objective Investigations}

  Based on the experiment results in previous section, hard instances for each algorithm are still hard for the others. In order to extract more useful knowledge about the strengths or weaknesses of certain algorithms on constraint algorithms, we need problem instances that are hard for one and easy for the others. Analysing the features of these instances helps us extracting knowledge regarding the strengths and weaknesses of algorithms by examining why an algorithm performs better on some groups of features while the others fails. This will help us developing more efficient prediction model for automated algorithm selection.

  To do this, we use a multi-objective DE algorithm (DEMO) described in \cite{robivc2005demo} to minimise the FEN for one algorithm and maximise it for the others. In other words, the FEN for generated problem instances is higher (harder) for a certain algorithm and lower (easier) for the others. In order to find instances that are hard for one algorithm type and easy for the others, we need to find solution as diverse as possible. Also, the solutions need to be close to pareto front. Satisfying these two aims makes us to use multi-objective evolutionary algorithm to generate problem instances. Hence, we use differential evolution for multi-objective optimisation (DEMO) proposed by Robic in \cite{robivc2005demo}. Based on results in \cite{robivc2005demo}, the DEMO achieves efficiently the above two goals. In DEMO, the candidate solution replaces parent when it dominates it and if the parent dominates it, the candidate is discarded. Otherwise, if the candidate and parent cannot dominate each other, the candidate is added to the population. The major difference between DEMO and other multi-objective evolutionary algorithms is that the newly generated good candidates are immediately used in creation of the subsequent candidates. This improves fast convergence to the true pareto front, while the use of non-dominated sorting and crowding distance metric in truncation of the extended population promotes the uniform
  spread of solutions. We refer the reader to \cite{robivc2005demo} for further details and implementation.

  In the following, we discuss the results for algorithms performances comparison. We carry out 30 independent runs for each number of constraints that are hard for one algorithm but still easy for the others. We set the evolving algorithm (DEMO) generation number to 5000 and the other parameters of evolving algorithm are set to pop size = $40$, CR = $0.5$, scaling factor = $0.9$ and $FEN_{max}$ is $300K$. Values for these parameters have been obtained by optimising the performance of the evolving algorithm in order to achieve the more easier and harder problem instances. For each of three algorithms, their best parameters are chosen \cite{takahama2010constrained,wang2009hybrid,hansen2010real}. First, the ($\epsilon$DEg) algorithm parameters are considered as: generation number = $1500$, pop size = $100$, CR = $0.5$, scaling factor = $0.5$. Also, the parameters for e-constraint method are described in \cite{shayan2015}. Moreover, for evolutionary strategy we perform $(1,7)$-ES algorithm with $1500$ generation using $P_{f} = 0.4$ with tendency to focus on feasible solution. In HMPSO algorithm, the swarm size $N$ is set to $60$, each sub-swarm size ($N_{s}$) is $8$ and all the PSO parameters are considered as Krohling and Coelho's PSO \cite{krohling2006coevolutionary}. In order to solve generated COPs, HMPSO generation number is set to $1500$. We need to say the parameters for the solvers are identical to those given in \cite{takahama2010constrained,wang2009hybrid,arnold20121+,robivc2005demo}

  In our all experiments, we generate set of problem instances that are hard to one algorithm and easy to the other ones. Tables \ref{table:fenmultilinsphere}, \ref{table:fenmultiquadsphere}, \ref{table:fenmultilinAckley}, \ref{table:fenmultiquadAckley}, \ref{table:fenmultilinRosen} and \ref{table:fenmultiquadRosen} show the function evaluation number (FEN) required for each algorithm to solve DE/ES/PSO hard instances for Sphere, Ackley and Rosenbrock objective functions (with 1 to 5 linear/quadratic constraints). As it is observed, there is more difference between the required FEN of instances generated by multi-objective algorithm evolver than the single-objective one.
  For instance, the required FEN for solving DE hard instances are higher for DE algorithm than solving it by ES and PSO algorithm. It means the DE hard instances are only hard for DE algorithm and easy for the others. In the following we start analysing constraint features of instances that are hard for one and easy for others.

\begin{table}
	\hspace{-1.0cm}
	\resizebox{0.45\columnwidth}{!}{%
		\parbox{.6\linewidth}{
			\centering

			\caption {The comparison of algorithms performance on each other's easy and hard instances based on required FEN for \textbf{Sphere} objective function with \textbf{linear} constraints. DE Easy (1 c) means instances that are easy for DE and with 1 constraint.}
			{%
				\begin{tabular}{|p{2.2cm}|c|c|c|c|}

					\hline \hline
					Instances	& DE algorithm & ES algorithm & PSO algorithm   \\ \hline
					DE Easy (1 c) &  25.6K &   28.2K & 33.2K  \\ \hline
					ES Easy (1 c) & 26.3K  & 27.1K   & 33.9K  \\ \hline
					PSO Easy (1 c) &  24.9K & 29.1K  & 72.5K   \\ \hline
					DE Easy (2 c) & 28.9K  & 21.9K   & 32.1K  \\ \hline
					ES Easy (2 c) &  25.2K &  24.3K  & 29.4K  \\ \hline
					PSO Easy (2 c) & 24.2K  & 25.2K   & 33.5K  \\ \hline
					DE Easy (3 c) & 32.4K  & 31.2K   & 33.9K  \\ \hline
					ES Easy (3 c) & 31.8K  &  29.1K  & 33.2K  \\ \hline
					PSO Easy (3 c) & 35.1K  & 28.6K   & 35.1K  \\ \hline
					DE Easy (4 c) & 34.2K  & 29.8K   & 38.2K  \\ \hline
					ES Easy (4 c) & 32.1K  & 31.5K   & 36.1K  \\ \hline
					PSO Easy (4 c) & 35.7K  & 28.9K   & 39.5K  \\ \hline
					DE Easy (5 c) &  35.3K &  42.1K  & 46.4K  \\ \hline
					ES Easy (5 c) &  31.2K &  45.2K  & 38.2K  \\ \hline
					PSO Easy (5 c) & 35.3K  & 44.9K   & 41.2K  \\ \hline \hline
					DE Hard (1 c) & 91.2K  & 78.3K   & 76.4K  \\ \hline
					ES Hard (1 c) &  81.3K &  86.4K  & 78.8K  \\ \hline
					PSO Hard (1 c) & 82.5K  & 72.5K   & 85.4K  \\ \hline
					DE Hard (2 c) & 93.4K &  81.3K  & 81.4K  \\ \hline
					ES Hard (2 c) & 84.3K  & 92.6K   & 79.4K  \\ \hline
					PSO Hard (2 c) &  85.7K & 84.1K   & 89.4K  \\ \hline
					DE Hard (3 c) & 98.3K  &  93.8K  & 78.9K  \\ \hline
					ES Hard (3 c) & 91.4K  & 108.6K  & 81.2K  \\ \hline
					PSO Hard (3 c) & 89.1K  & 98.2K   & 91.6K  \\ \hline
					DE Hard (4 c) & 104.2K  & 89.4K   & 82.5K  \\ \hline
					ES Hard (4 c) & 89.4K  & 115.1K   &  78.4K \\ \hline
					PSO Hard (4 c) &  92.9K &  93.5K  &  115.3K \\ \hline
					DE Hard (5 c) &  123.2K & 111.4K   & 98.4K  \\ \hline
					ES Hard (5 c) & 98.2K  &  133.2K  & 94.9K  \\ \hline
					PSO Hard (5 c) &  101.3K & 109.2K & 118.3K   \\ \hline
				\end{tabular}}
				\setlength{\tabcolsep}{5em}

				\label{table:perflinearsphere}

			}
		}
		\hfill
		\resizebox{0.45\columnwidth}{!}{%
			\parbox{.6\linewidth}{
				\centering
				\caption {The comparison of algorithms performance on each other's easy and hard instances based on required FEN for \textbf{Sphere} objective function with \textbf{quadratic} constraints. DE Easy (1 c) means instances that are easy for DE and with 1 constraint.}
				{%
					\begin{tabular}{|p{2.2cm}|c|c|c|c|}
					\hline \hline
					Instances	& DE algorithm & ES algorithm & PSO algorithm   \\ \hline
					DE Easy (1 c) & 24.2K  & 23.6K   & 24.9K  \\ \hline
					ES Easy (1 c) &  24.8K &  24.2K  & 25.4K  \\ \hline
					PSO Easy (1 c) & 26.4K  &  25.4K  & 26.4K  \\ \hline
					DE Easy (2 c) & 25.3K  & 28.1K   & 26.4K  \\ \hline
					ES Easy (2 c) & 24.1K  &  27.2K  & 27.4K  \\ \hline
					PSO Easy (2 c) & 23.5K  &  29.3K  & 271.K  \\ \hline
					DE Easy (3 c) &  27.9K & 31.9K   & 35.5K  \\ \hline
					ES Easy (3 c) &  29.4K &  32.1K  &  28.5K \\ \hline
					PSO Easy (3 c) & 28.1K  & 28.7K   &  29.4K \\ \hline
					DE Easy (4 c) & 34.1K  & 28.9K   & 36.4K  \\ \hline
					ES Easy (4 c) & 35.2K  & 35.3K   &  31.6K \\ \hline
					PSO Easy (4 c) &  31.8K &  29.5K  & 33.2K  \\ \hline
					DE Easy (5 c) &  38.7K & 29.2K   &  37.2K \\ \hline
					ES Easy (5 c) &  35.6K &  28.2K  & 39.5K  \\ \hline
					PSO Easy (5 c) & 36.3K  & 31.5K   & 36.2K  \\ \hline \hline
					DE Hard (1 c) & 129.3K  & 102.7K   & 105.3K  \\ \hline
					ES Hard (1 c) &   104.3K & 121.2K   &  108.2K \\ \hline
					PSO Hard (1 c) & 108.2K  & 104.2K   & 119.8K  \\ \hline
					DE Hard (2 c) &  132.6K & 114.2K   & 114.9K  \\ \hline
					ES Hard (2 c) & 111.2K  & 127.1K   & 112.4K  \\ \hline
					PSO Hard (2 c) & 109.4K  & 112.4K   & 125.3K  \\ \hline
					DE Hard (3 c) & 136.2K  & 116.3K   & 112.4K  \\ \hline
					ES Hard (3 c) & 117.2K  &  132.1K  &  109.9K \\ \hline
					PSO Hard (3 c) & 119.8K  & 119.2K   & 132.6K  \\ \hline
					DE Hard (4 c) & 141.2K  & 119.9K   & 119.6K  \\ \hline
					ES Hard (4 c) & 113.8K  & 131.2K   &  121.9K \\ \hline
					PSO Hard (4 c) & 115.4K  & 121.4K   & 138.9K  \\ \hline
					DE Hard (5 c) & 149.3K  & 129.7K   &  122.9K \\ \hline
					ES Hard (5 c) &  124.4K & 149.6K   &  126.4K \\ \hline
					PSO Hard (5 c) &  123.9K &  124.2K  & 148.3K  \\ \hline
					\end{tabular}}
				\setlength{\tabcolsep}{5em}

					\label{table:perfquadsphere}
				}
			}

		\end{table}

	\begin{table}
		\hspace{-1.0cm}
		\resizebox{0.45\columnwidth}{!}{%
			\parbox{.6\linewidth}{
				\centering

				\caption {The comparison of algorithms performance on each other's easy and hard instances based on required FEN for \textbf{Achkley} objective function with \textbf{linear} constraints. DE Easy (1 c) means instances that are easy for DE and with 1 constraint.}
				{%
					\begin{tabular}{|p{2.2cm}|c|c|c|c|}

						\hline \hline
						Instances	& DE algorithm & ES algorithm & PSO algorithm   \\ \hline
						DE Easy (1 c) &  39.2K & 37.1K  & 37.4K   \\ \hline
						ES Easy (1 c) &  38.6K &  37.8K  & 38.1K  \\ \hline
						PSO Easy (1 c) & 41.3K  & 39.2K   & 41.7K   \\ \hline
						DE Easy (2 c) & 40.3K  &  41.5K  & 42.6K  \\ \hline
						ES Easy (2 c) &  41.3K & 42.5K   & 41.5K  \\ \hline
						PSO Easy (2 c) & 42.6K  &  41.2K  & 44.7K  \\ \hline
						DE Easy (3 c) &  47.3K & 48.2K   & 46.9K  \\ \hline
						ES Easy (3 c) &  48.9K &  47.3K  & 46.3K  \\ \hline
						PSO Easy (3 c) & 49.2K  &  51.6K  & 50.9K  \\ \hline
						DE Easy (4 c) & 48.9K  & 47.2K   & 49.2K  \\ \hline
						ES Easy (4 c) &  49.2K &  48.2K  & 50.1K  \\ \hline
						PSO Easy (4 c) & 51.2K  & 50.5K   & 52.6K  \\ \hline
						DE Easy (5 c) & 51.3K  & 52.7K   & 51.6K  \\ \hline
						ES Easy (5 c) &  52.1K & 52.6K   & 50.7K  \\ \hline
						PSO Easy (5 c) & 55.3K  & 54.7K   & 52.3K  \\ \hline
						\hline
						DE Hard (1 c) & 107.3K  &  83.2K  & 85.6K  \\ \hline
						ES Hard (1 c) &  82.3K &  105.2K  & 81.6K  \\ \hline
						PSO Hard (1 c) &  85.2K &  83.9K  & 106.3K  \\ \hline
						DE Hard (2 c) & 114.2K  & 88.2K   & 91.5K  \\ \hline
						ES Hard (2 c) &  87.3K &  115.3K  & 88.8K  \\ \hline
						PSO Hard (2 c) &  89.2K &   87.3K &  116.9K \\ \hline
						DE Hard (3 c) &  119.8K &  94.1K  & 93.9K  \\ \hline
						ES Hard (3 c) &  95.2K &  121.6K  &  94.2K \\ \hline
						PSO Hard (3 c) &  93.2K & 95.1K   & 121.5K  \\ \hline
						DE Hard (4 c) & 125.2K  & 99.2K   & 101.4K  \\ \hline
						ES Hard (4 c) &  101.3K &  126.3K  & 98.2K  \\ \hline
						PSO Hard (4 c) &  99.4K &  97.8K  &  127.4K \\ \hline
						DE Hard (5 c) & 132.5K  & 102.2K   & 101.5K  \\ \hline
						ES Hard (5 c) &  101.4K &  134.7K  & 103.5K  \\ \hline
						PSO Hard (5 c) &  103.9K &  102.4K  & 131.4K  \\ \hline
					\end{tabular}}
					\setlength{\tabcolsep}{5em}

					\label{table:perflinearackley}

				}
			}
			\hfill
			\resizebox{0.45\columnwidth}{!}{%
				\parbox{.6\linewidth}{
					\centering
					\caption {The comparison of algorithms performance on each other's easy and hard instances based on required FEN for \textbf{Ackley} objective function with \textbf{quadratic} constraints. DE Easy (1 c) means instances that are easy for DE and with 1 constraint.}
					{%
						\begin{tabular}{|p{2.2cm}|c|c|c|c|}
							\hline \hline
							Instances	& DE algorithm & ES algorithm & PSO algorithm   \\ \hline
							DE Easy (1 c) &  38.1K &  39.4K  & 37.2K  \\ \hline
							ES Easy (1 c) &  37.1K & 38.1K   & 39.0K  \\ \hline
							PSO Easy (1 c) &  41.2K & 39.9K   & 40.7K  \\ \hline
							DE Easy (2 c) & 38.9K  &  43.9K  & 41.7K  \\ \hline
							ES Easy (2 c) & 40.1K  & 41.2K   & 43.2K  \\ \hline
							PSO Easy (2 c) &  39.1K &   43.1K &  46.1K \\ \hline
							DE Easy (3 c) & 46.3K  & 47.9K   & 45.1K  \\ \hline
							ES Easy (3 c) & 49.1K  &  48.1K  & 47.2K  \\ \hline
							PSO Easy (3 c) &  42.1K & 45.1K   & 49.1K  \\ \hline
							DE Easy (4 c) &  49.2K &  48.7K  &  48.4K \\ \hline
							ES Easy (4 c) &  49.7K &  49.9K  &   52.9K\\ \hline
							PSO Easy (4 c) & 56.1K  & 55.1K   & 54.1K  \\ \hline
							DE Easy (5 c) &  50.0K & 51.2K   & 52.4K  \\ \hline
							ES Easy (5 c) &  51.3K & 56.2K   & 54.1K  \\ \hline
							PSO Easy (5 c) & 61.2K  & 58.9K   & 59.1K  \\ \hline \hline
							DE Hard (1 c) & 133.4K  &  93.1K  & 94.6K  \\ \hline
							ES Hard (1 c) & 92.1K  &  135.1K  & 94.1K  \\ \hline
							PSO Hard (1 c) & 95.2K  & 94.9K   & 134.2K  \\ \hline
							DE Hard (2 c) & 139.1K  & 98.2K   & 97.1K  \\ \hline
							ES Hard (2 c) &  97.1K &   138.2K &  99.1K \\ \hline
							PSO Hard (2 c) &  109.1K &  107.2K  & 145.2K  \\ \hline
							DE Hard (3 c) &  145.3K &  112.6K  &  109.6K \\ \hline
							ES Hard (3 c) &  111.6K  &  141.2K  & 108.9K  \\ \hline
							PSO Hard (3 c) &  111.2K &  109.2K  & 152.K  \\ \hline
							DE Hard (4 c) & 167.2K  & 132.1K   & 135.1K  \\ \hline
							ES Hard (4 c) & 131.1K   & 167.9K   &  133.9K \\ \hline
							PSO Hard (4 c) & 132.1K  & 133.2K   & 169.1K  \\ \hline
							DE Hard (5 c) & 177.1K   & 143.2K   & 131.3K  \\ \hline
							ES Hard (5 c) &  141.2K &  181.2K  & 144.9K  \\ \hline
							PSO Hard (5 c) & 139.1K  &  142.9K  & 182.1K  \\ \hline
						\end{tabular}}
						\setlength{\tabcolsep}{5em}

						\label{table:perfquadackley}
					}
				}

			\end{table}
\begin{table}
	\hspace{-1.0cm}
	\resizebox{0.45\columnwidth}{!}{%
		\parbox{.6\linewidth}{
			\centering

			\caption {The comparison of algorithms performance on each other's easy and hard instances based on required FEN for \textbf{Rosenbrock} objective function with \textbf{linear} constraints. DE Easy (1 c) means instances that are easy for DE and with 1 constraint.}
			{%
				\begin{tabular}{|p{2.2cm}|c|c|c|c|}

					\hline \hline
					Instances	& DE algorithm & ES algorithm & PSO algorithm   \\ \hline
					DE Easy (1 c) &  38.1K & 37.4K   & 39.1K  \\ \hline
					ES Easy (1 c) &  39.2K & 38.9K   & 36.2K  \\ \hline
					PSO Easy (1 c) & 44.6K  &   40.1K & 43.1K  \\ \hline
					DE Easy (2 c) & 41.2K  & 42.4K   & 47.1K  \\ \hline
					ES Easy (2 c) & 40.2K  &  45.2K  & 40.4K  \\ \hline
					PSO Easy (2 c) &  43.9K &   42.6K & 44.8K  \\ \hline
					DE Easy (3 c) &  41.2K & 42.3K   & 45.4K  \\ \hline
					ES Easy (3 c) &  47.7K &   48.1K & 47.4K  \\ \hline
					PSO Easy (3 c) & 48.3K  &  52.4K  & 53.1K  \\ \hline
					DE Easy (4 c) & 44.1K  & 42.7K   & 43.3K  \\ \hline
					ES Easy (4 c) &  48.4K &  49.7K  & 52.6K  \\ \hline
					PSO Easy (4 c) & 53.5K  &  53.1K  &  54.9K \\ \hline
					DE Easy (5 c) & 55.5K  & 53.2K   &  52.1K \\ \hline
					ES Easy (5 c) &  52.8K & 55.4K   &  52.1K \\ \hline
					PSO Easy (5 c) & 58.2K  &   53.4K &  52.8K \\ \hline \hline
					DE Hard (1 c) &  110.2K &  83.7K  & 85.2K  \\ \hline
					ES Hard (1 c) &  81.5K &  107.2K  & 82.1K  \\ \hline
					PSO Hard (1 c) &  86.9K &  85.3K  & 108.2K  \\ \hline
					DE Hard (2 c) & 116.6K  & 89.3K   & 92.1K  \\ \hline
					ES Hard (2 c) &  88.2K &  117.5K  & 87.0K  \\ \hline
					PSO Hard (2 c) &  90.8K &  88.1K  & 114.5K  \\ \hline
					DE Hard (3 c) & 121.2K  & 92.1K   & 92.5K  \\ \hline
					ES Hard (3 c) & 94.1K  &  124.8K  &  93.2K \\ \hline
					PSO Hard (3 c) & 94.7K  &  96.5K  &  125.2K \\ \hline
					DE Hard (4 c) &  126.1K &  101.6K  & 104.2K  \\ \hline
					ES Hard (4 c) & 104.8K  & 127.2K   & 96.1K  \\ \hline
					PSO Hard (4 c) &  100.2K &  92.1K  & 123.7K  \\ \hline
					DE Hard (5 c) &  135.1K &  109.5K  &  106.8K \\ \hline
					ES Hard (5 c) &  105.2K &   136.1K &   105.1K \\ \hline
					PSO Hard (5 c) &  102.1K & 106.8K   & 131.4  \\ \hline
				\end{tabular}}
				\setlength{\tabcolsep}{5em}

				\label{table:perflinearrosenbrock}

			}
		}
		\hfill
		\resizebox{0.45\columnwidth}{!}{%
			\parbox{.6\linewidth}{
				\centering
				\caption {The comparison of algorithms performance on each other's easy and hard instances based on required FEN for \textbf{Rosenbrock} objective function with \textbf{quadratic} constraints. DE Easy (1 c) means instances that are easy for DE and with 1 constraint.}
				{%
					\begin{tabular}{|p{2.2cm}|c|c|c|c|}
						\hline \hline
						Instances	& DE algorithm & ES algorithm & PSO algorithm   \\ \hline
						DE Easy (1 c) & 41.2K  & 40.2K   & 38.7K  \\ \hline
						ES Easy (1 c) &  39.2K & 39.3K   & 36.1K  \\ \hline
						PSO Easy (1 c) &  43.1K & 39.6K   & 42.2K  \\ \hline
						DE Easy (2 c) & 39.1K  & 44.4K   & 43.2K  \\ \hline
						ES Easy (2 c) &   42.6K & 44.5K   & 41.8K  \\ \hline
						PSO Easy (2 c) & 41.2K & 45.6K  & 47.2K    \\ \hline
						DE Easy (3 c) &  47.1K &  48.5K  &  49.2K \\ \hline
						ES Easy (3 c) &  46.8K & 49.5K   & 48.8K   \\ \hline
						PSO Easy (3 c) &  46.2K &  42.7K  & 48.4K  \\ \hline
						DE Easy (4 c) &  48.7K &  49.1K  & 51.2K  \\ \hline
						ES Easy (4 c) &  50.2K &  52.4K  & 55.2K  \\ \hline
						PSO Easy (4 c) &  59.2K &  54.5K  &  51.9K \\ \hline
						DE Easy (5 c) &  52.5K &  56.1K  & 55.7K  \\ \hline
						ES Easy (5 c) &  56.0K &  55.7K  & 53.8K  \\ \hline
						PSO Easy (5 c) & 66.3K  & 59.8K   & 60.4K  \\ \hline \hline
						DE Hard (1 c) & 137.2K  & 93.7K   & 92.1K   \\ \hline
						ES Hard (1 c) & 93.7K  & 138.2K   & 99.2K  \\ \hline
						PSO Hard (1 c) & 93.1K  &  96.2K  &  138.0K \\ \hline
						DE Hard (2 c) & 142.7K  & 99.7K   & 98.4K  \\ \hline
						ES Hard (2 c) &  98.8K &  141.6K  & 101.4K  \\ \hline
						PSO Hard (2 c) &  112.7K &  109.5K  & 148.1K  \\ \hline
						DE Hard (3 c) & 148.2K  & 115.3K   & 112.8K  \\ \hline
						ES Hard (3 c) &  114.1K &  144.8K  &  113.2K \\ \hline
						PSO Hard (3 c) & 113.6K  & 113.1K   & 157.3K   \\ \hline
						DE Hard (4 c) &  171.7K &  136.7K  & 133.4K  \\ \hline
						ES Hard (4 c) &  134.7K &  168.2K  &  135.3K \\ \hline
						PSO Hard (4 c) &  134.7K & 139.3K   & 172.6K  \\ \hline
						DE Hard (5 c) &  179.6K &  146.1K  & 144.8K  \\ \hline
						ES Hard (5 c) &  143.8K &   185.1K &  147.4K \\ \hline
						PSO Hard (5 c) & 143.7K  & 141.4K   & 186.4K  \\ \hline
					\end{tabular}}
					\setlength{\tabcolsep}{5em}

					\label{table:perfquadrosenbrock}
				}
			}

		\end{table}

\begin{table}
	\hspace{-1.0cm}
	\resizebox{0.45\columnwidth}{!}{%
		\parbox{.6\linewidth}{
			\centering

			\caption {The FEN required for each algorithm to solve  DE/ES/PSO hard instances (Sphere for 1 to 5 \textbf{linear constraints})}
			{%
				\begin{tabular}{|p{2.1cm}|c|c|c|c|}

					\hline \hline
					Instances	& DE algorithm & ES algorithm & PSO algorithm   \\ \hline
					DE hard (1 c) & \textbf{86.3K}  &  41.5K  &  43.2K \\ \hline
					ES hard (1 c) & 45.7K & \textbf{84.2K}  & 48.3K    \\ \hline
					PSO hard (1 c)& 37.2K &  41.8K &  \textbf{80.1K}   \\ \hline
					DE hard (2 c) & \textbf{88.8K}  &  43.9K  &  44.2K\\ \hline
					ES hard (2 c) & 45.9K & \textbf{85.4K}  & 46.3K   \\ \hline
					PSO hard (2 c)& 43.2K &  42.5K &  \textbf{82.9K} \\	\hline
					DE hard (3 c) & \textbf{91.4K}  &  44.6K  &  45.3K \\ \hline
					ES hard (3 c) & 49.2K & \textbf{87.8K}  & 48.1K    \\ \hline
					PSO hard (3 c)& 46.2K &  47.7K &  \textbf{85.5K}   \\ \hline
					DE hard (4 c) & \textbf{94.2K}  &  47.5K  &  47.8K\\ \hline
					ES hard (4 c) & 51.7K & \textbf{89.1K}  & 50.1K   \\ \hline
					PSO hard (4 c)& 48.7K &  49.9K &  \textbf{87.3K} \\ \hline
					DE hard (5 c) & \textbf{96.2K}  &  48.2K  &  49.5K \\ \hline
					ES hard (5 c) & 52.4K & \textbf{90.4K}  & 53.5K    \\ \hline
					PSO hard (5 c)& 49.6K &  51.4K &  \textbf{91.6K}   \\ \hline
				\end{tabular}}
				\setlength{\tabcolsep}{5em}

				\label{table:fenmultilinsphere}

			}
		}
		\hfill
		\resizebox{0.45\columnwidth}{!}{%
			\parbox{.6\linewidth}{
				\centering
				\caption {TThe FEN required for each algorithm to solve  DE/ES/PSO hard instances (Sphere for 1 to 5 \textbf{quadratic constraints})}
				{%
					\begin{tabular}{|p{2.1cm}|c|c|c|c|}

						\hline \hline
						Instances	& DE algorithm & ES algorithm & PSO algorithm   \\ \hline
						DE hard (1 c) & \textbf{92.3K}  &  50.2K  &  51.9K \\ \hline
						ES hard (1 c) & 48.8K & \textbf{91.3K}  & 49.3K    \\ \hline
						PSO hard (1 c)& 44.5K &  46.8K &  \textbf{93.1K}   \\ \hline
						DE hard (2 c) & \textbf{93.5K}  &  52.9K  &  54.2K\\ \hline
						ES hard (2 c) & 50.9K & \textbf{95.9K}  & 51.2K   \\ \hline
						PSO hard (2 c)& 50.2K &  53.2K &  \textbf{96.3K} \\	\hline
						DE hard (3 c) & \textbf{95.9K}  &  54.3K  &  55.3K \\ \hline
						ES hard (3 c) & 53.9K & \textbf{97.4K}  & 52.4K    \\ \hline
						PSO hard (3 c)& 57.3K &  56.3K &  \textbf{98.9K}   \\ \hline
						DE hard (4 c) & \textbf{98.3K}  &  56.4K  &  57.3K\\ \hline
						ES hard (4 c) & 56.3K & \textbf{102.3K}  & 52.1K   \\ \hline
						PSO hard (4 c)& 59.2K &  58.2K &  \textbf{101.6K} \\ \hline
						DE hard (5 c) & \textbf{102.1K}  &  58.3K  &  59.4K \\ \hline
						ES hard (5 c) & 59.2K & \textbf{103.2K}  & 60.2K    \\ \hline
						PSO hard (5 c)& 62.6K &  63.8K &  \textbf{105.2K}   \\ \hline
					\end{tabular}}
					\setlength{\tabcolsep}{5em}

					\label{table:fenmultiquadsphere}
				}
			}

		\end{table}

\begin{table}
	\hspace{-1.0cm}
	\resizebox{0.45\columnwidth}{!}{%
		\parbox{.6\linewidth}{
			\centering

			\caption {The FEN required for each algorithm to solve  DE/ES/PSO hard instances (Ackley for 1 to 5 \textbf{linear constraints})}
			{%
				\begin{tabular}{|p{2.1cm}|c|c|c|c|}

					\hline \hline
					Instances	& DE algorithm & ES algorithm & PSO algorithm   \\ \hline
					DE hard (1 c) & \textbf{102.3K}  &  46.1K  &  51.4K \\ \hline
					ES hard (1 c) & 51.2K & \textbf{104.7K}  & 50.2K    \\ \hline
					PSO hard (1 c)& 47.4K &  49.8K &  \textbf{107.4K}   \\ \hline
					DE hard (2 c) & \textbf{112.1K}  &  56.1K  &  54.1K\\ \hline
					ES hard (2 c) & 53.9K & \textbf{115.9K}  & 48.6K   \\ \hline
					PSO hard (2 c)& 55.5K &  55.3K &  \textbf{117.2K} \\	\hline
					DE hard (3 c) & \textbf{126.1K}  &  63.7K  &  65.2K \\ \hline
					ES hard (3 c) & 59.1K & \textbf{128.3K}  & 58.7K    \\ \hline
					PSO hard (3 c)& 61.7K &  62.8K &  \textbf{134.2K}   \\ \hline
					DE hard (4 c) & \textbf{124.9K}  &  68.4K  &  63.1K\\ \hline
					ES hard (4 c) & 64.1K & \textbf{129.8K}  & 59.2K   \\ \hline
					PSO hard (4 c)& 67.5K &  69.2K &  \textbf{135.2K} \\ \hline
					DE hard (5 c) & \textbf{138.8K}  &  75.2K  &  74.1K \\ \hline
					ES hard (5 c) & 71.2K & \textbf{137.1K}  & 76.7K    \\ \hline
					PSO hard (5 c)& 73.1K &  74.1K &  \textbf{141.2K}   \\ \hline
				\end{tabular}}
				\setlength{\tabcolsep}{5em}

				\label{table:fenmultilinAckley}

			}
		}
		\hfill
		\resizebox{0.45\columnwidth}{!}{%
			\parbox{.6\linewidth}{
				\centering
				\caption {The FEN required for each algorithm to solve  DE/ES/PSO hard instances (Ackley for 1 to 5 \textbf{quadratic constraints})}
				{%
					\begin{tabular}{|p{2.1cm}|c|c|c|c|}

						\hline \hline
						Instances	& DE algorithm & ES algorithm & PSO algorithm   \\ \hline
						DE hard (1 c) & \textbf{142.5K}  &  60.1K  &  62.5K \\ \hline
						ES hard (1 c) & 58.5K & \textbf{148.2K}  & 61.4K    \\ \hline
						PSO hard (1 c)& 53.2K &  53.9K &  \textbf{147.7K}   \\ \hline
						DE hard (2 c) & \textbf{153.3K}  &  58.1K  &  58.1K\\ \hline
						ES hard (2 c) & 59.2K & \textbf{155.5K}  & 59.2K   \\ \hline
						PSO hard (2 c)& 57.8K &  56.3K &  \textbf{157.2K} \\	\hline
						DE hard (3 c) & \textbf{167.3K}  &  65.2K  &  68.1K \\ \hline
						ES hard (3 c) & 63.2K & \textbf{169.2K}  & 69.8K    \\ \hline
						PSO hard (3 c)& 65.7K &  67.9K &  \textbf{167.6K}   \\ \hline
						DE hard (4 c) & \textbf{174.8K}  &  71.2K  &  75.1K\\ \hline
						ES hard (4 c) & 66.8K & \textbf{169.1K}  & 72.9K   \\ \hline
						PSO hard (4 c)& 69.1K &  68.3K &  \textbf{172.9K} \\ \hline
						DE hard (5 c) & \textbf{179.5K}  &  75.1K  &  76.1K \\ \hline
						ES hard (5 c) & 72.8K & \textbf{174.9K}  & 77.4.2K    \\ \hline
						PSO hard (5 c)& 75.1K &  74.9K &  \textbf{175.9K}   \\ \hline
					\end{tabular}}
					\setlength{\tabcolsep}{5em}

					\label{table:fenmultiquadAckley}
				}
			}

		\end{table}

\begin{table}
	\hspace{-1.0cm}
	\resizebox{0.45\columnwidth}{!}{%
		\parbox{.6\linewidth}{
			\centering

			\caption {The FEN required for each algorithm to solve  DE/ES/PSO hard instances (Rosenbrock for 1 to 5 \textbf{linear constraints})}
			{%
				\begin{tabular}{|p{2.1cm}|c|c|c|c|}

					\hline \hline
					Instances	& DE algorithm & ES algorithm & PSO algorithm   \\ \hline
					DE hard (1 c) & \textbf{103.1K}  &  48.2K  &  53.9K \\ \hline
					ES hard (1 c) & 53.1K & \textbf{107.2K}  & 52.7K    \\ \hline
					PSO hard (1 c)& 45.7K &  48.1K &  \textbf{109.2K}   \\ \hline
					DE hard (2 c) & \textbf{115.1K}  &  57.8K  &  55.7K\\ \hline
					ES hard (2 c) & 54.7K & \textbf{113.4K}  & 46.1K   \\ \hline
					PSO hard (2 c)& 54.8K &  54.9K &  \textbf{119.5K} \\	\hline
					DE hard (3 c) & \textbf{124.3K}  &  65.2K  &  62.1K \\ \hline
					ES hard (3 c) & 58.8K & \textbf{127.1K}  & 59.7K    \\ \hline
					PSO hard (3 c)& 62.3K &  65.5K &  \textbf{136.1K}   \\ \hline
					DE hard (4 c) & \textbf{125.5K}  &  69.1K  &  65.2K\\ \hline
					ES hard (4 c) & 67.5K & \textbf{128.1K}  & 58.7K   \\ \hline
					PSO hard (4 c)& 64.9K &  70.6K &  \textbf{137.1K} \\ \hline
					DE hard (5 c) & \textbf{135.1K}  &  74.1K  &  74.7K \\ \hline
					ES hard (5 c) & 73.7K & \textbf{135.8K}  & 75.1K    \\ \hline
					PSO hard (5 c)& 72.3K &  76.5K &  \textbf{140.9K}   \\ \hline
				\end{tabular}}
				\setlength{\tabcolsep}{5em}

				\label{table:fenmultilinRosen}

			}
		}
		\hfill
		\resizebox{0.45\columnwidth}{!}{%
			\parbox{.6\linewidth}{
				\centering
				\caption {The FEN required for each algorithm to solve  DE/ES/PSO hard instances (Rosenbrock for 1 to 5 \textbf{quadratic constraints})}
				{%
					\begin{tabular}{|p{2.1cm}|c|c|c|c|}

						\hline \hline
						Instances	& DE algorithm & ES algorithm & PSO algorithm   \\ \hline
						DE hard (1 c) & \textbf{143.1K}  &  61.4K  &  63.7K \\ \hline
						ES hard (1 c) & 59.6K & \textbf{149.7K}  & 62.5K    \\ \hline
						PSO hard (1 c)& 54.3K &  54.2K &  \textbf{143.8K}   \\ \hline
						DE hard (2 c) & \textbf{155.2K}  &  59.2K  &  59.2K\\ \hline
						ES hard (2 c) & 61.3K & \textbf{154.8K}  & 57.9K   \\ \hline
						PSO hard (2 c)& 59.2K &  57.1K &  \textbf{158.0K} \\	\hline
						DE hard (3 c) & \textbf{168.9K}  &  63.7K  &  66.8K \\ \hline
						ES hard (3 c) & 65.2K & \textbf{170.1K}  & 68.1K    \\ \hline
						PSO hard (3 c)& 63.7K &  68.1K &  \textbf{168.9K}   \\ \hline
						DE hard (4 c) & \textbf{175.1K}  &  73.8K  &  76.9K\\ \hline
						ES hard (4 c) & 68.2K & \textbf{172.7K}  & 75.2K   \\ \hline
						PSO hard (4 c)& 67.7K &  69.1K &  \textbf{176.2K} \\ \hline
						DE hard (5 c) & \textbf{180.2K}  &  74.2K  &  77.8K \\ \hline
						ES hard (5 c) & 74.2K & \textbf{175.1K}  & 79.2K    \\ \hline
						PSO hard (5 c)& 73.6K &  74.4K &  \textbf{179.4K}   \\ \hline
					\end{tabular}}
					\setlength{\tabcolsep}{5em}

					\label{table:fenmultiquadRosen}
				}
			}

		\end{table}
\subsection{Analysis for Linear Constraints }

We run our experiments on Sphere, Ackley and Rosenbrock objective functions. The linear constraints are considered as in Equation \ref{lineareq} with all coefficients $a_{n}$s that are in the range of $[-5,5]$. Also, the problem dimension is set to 30. As it mentioned before, to analyse and discuss some features such as shortest distance, we assume that the optimum is zero ($b\leq 0$). We use three types of problem instances. DE hard denotes problem instances that are hard for DE algorithm but still easy for PSO and ES algorithms. Also, ES hard instances are easy for DE and PSO algorithm in this section. PSO hard means the instances that are hard for PSO but easy for the rest. Each constraint is generated using multi-objective evolver to generate instances that are hard for one algorithm but easy for others. In the following we discuss the features of linear constraints.

Figure \ref{fig:box_plot_std} represents some evidence of linear constraint coefficient relationship (standard deviation). It is shown that standard deviation of (1 to 5) linear constraints are higher for DE hard instances than ES and PSO hard ones. This result is similar for all Sphere, Ackley and Rosenbrock objective functions. This means, the instances that are hard for DE algorithm but easy for ES and PSO have higher standard deviation for their constraints coefficients. In other words, this constraint feature has influence on problem difficulty. This improves the prediction ability for algorithm selection framework.

Box plots shown in Figure \ref{fig:box_plot_dis} represent the shortest distance from optimum feature for hard instances. Based on the experiments, hard instances for ES algorithm have higher value (closer to optimum) shortest distance than the other algorithms. It is noteworthy that lower value in Figure \ref{fig:box_plot_dis} means the constraint hyperplane is further from optimum. In other words, the constraints hyperplanes are closer to the optimum in ES hard instances. This relationship holds the pattern for all objective functions in linear constraints.
We also study the feasibility ratio in vicinity of the optimum. As observed in Table \ref{table:feasibilityratiolin}, hard DE instances have lower feasibility ratio comparing to PSO and ES hard instances. This follows the same pattern for all experimented objective functions. Also increasing the number of constraints decreases the problem optimum-local feasibility for all algorithm problem instances. The angle between linear constraints feature is analysed for linear constraints. As it is observed in Table \ref{table:anglelinear}, ES hard instances have lower angle values for all Sphere, Ackley and Rosenbrock objective functions. This means, instances that are hard for ES have less angle value between their constraint hyperplanes. Interestingly, all objective function that we use in this experiment follow the same relationship.

As it is observed, to compare the instances, DE hard instances have higher linear constraint coefficient standard deviation. It can be translated as DE algorithm has more difficulty to coefficients standard deviation feature than PSO and ES algorithms. Also, the local-optimum feasibility ratio value is higher in ES and PSO hard instances than DE hard ones. This means, ES and PSO algorithms are more effective to problems with higher optimum feasibility ratio feature. The shortest distance and angle features for ES is less than DE and PSO hard instances. Interestingly, this features are similar for all used objective functions. The linear constraint feature based analysis gives us helpful knowledge to implement algorithm selection framework.

\subsection{Analysis for Quadratic Constraints }
In this section, we carry out our experiments on Sphere, Ackley and Rosenbrock objective function with quadratic constraints (see Equation \ref{quadraticeq}) using same setup as previous section. In the following we do feature based analysis of constraints in hard DE, PSO and ES instances (that are easy for the other algorithms).

Figure \ref{fig:box_plot_std} shows some evidence of quadratic constraint coefficients relationship. Based on our experiments, in each constraint, the quadratic coefficient has more ability than linear coefficients to make problem harder to solve. In other words, in Equation \ref{quadraticeq}, $a_{1}$ is more contributing than $a_{2}$ to problem difficulty. As it is shown in the box plots, the standard deviation of 1 to 5 quadratic constraints in DE hard instances are higher comparing the other two algorithm hard instances. In contrast, our results show no systematic relationship between problem difficulty and linear coefficients in each quadratic constraints and quadratic coefficients have more contribution in problem difficulty.

As it is observed in Figure \ref{fig:box_plot_dis}, the shortest distance feature for DE, PSO and ES hard instances are compared. In instances that are hard for ES and easy for the other algorithms, the quadratic constraint hyperplanes are closer to optimum (zero). This applies to all experimented objective functions.
Also, calculating the angle feature for quadratic constraint does not show any systematic relationship to problem difficulty. The feasibility ratio near the optimum is analysed for DE, ES and PSO hard instances. As it is shown in Table \ref{table:feasibilityratioquad}, the feasibility ratio in DE hard instances are lower than the other algorithms hard instances. All objective functions have the same pattern. Also, the number of constraint has a systematic relationship with feasibility ratio.

Based on the results, to compare COP instances with quadratic constraints, DE hard instances have higher coefficient standard deviation value than the other algorithm hard ones. It is translated as the DE algorithm has more difficulty solving instances with higher standard deviation value for their quadratic constraints than ES and PSO. Also, the quadratic constraints are closer to optimum in ES instances than the other experimented algorithms. In other words, ES algorithm is more influenced by constraint with closer to optimum instances. Moreover, the optimum feasibility ratio in DE instances are lower than PSO and ES.\\

\begin{table}[t]
	\centering
	\caption {The angle feature for Sphere objective function for linear constraints}
	\resizebox{\columnwidth}{1.3 cm}{%
		\begin{tabular}{|p{1.5cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|c|}

			\hline \hline
			& Cons 1,2 & Cons 1,3 & Cons 1,4 & Cons 1,5 & Cons 2,3 & Cons 2,4 & Cons 2,5 & Cons 3,4 & Cons 3,5 & Cons 4,5 \\ \hline

			DE Hard & 74  & 64 &   63  & 58  & 74 & 71  & 68  &   59  & 62  & 86\\ \hline
			ES Hard 	& 33  & 21  &  37  & 24  & 44 & 46  & 39  &   46  & 48  & 51\\ \hline
			PSO Hard & 75  & 63  &   82 & 68  & 71 & 73  & 72  &   69  & 81  & 86 \\ \hline

		\end{tabular}}
		\setlength{\tabcolsep}{5em}

		\label{table:anglelinear}
	\end{table}

	\begin{table}[t]
		\caption {Optimum-local feasibility ratio of search space near the optimum for 1,2,3,4 and 5 linear constraint}
		\centering
		\resizebox{0.65\columnwidth}{!}{%

			\begin{tabular}{|p{1.5cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|c|}

				\hline \hline
				& 1 cons & 2 cons & 3 cons & 4 cons & 5 cons   \\ \hline
				DE Hard & 6\% & 5\%  &  3\%   &  3 \%  & 2\%  \\ \hline
				ES Hard & 16\% &  11\% &    10 \%  &   6\% &  5\% \\ \hline
				PSO Hard & 17\% & 12\%  &   11\%  &   8\% &  5\%  \\ \hline

			\end{tabular}}
			\setlength{\tabcolsep}{5em}

			\label{table:feasibilityratiolin}
		\end{table}

		\begin{table}[t]
			\caption {Optimum-local feasibility ratio of search space near the optimum for 1,2,3,4 and 5 quadratic constraint}
			\centering
			\resizebox{0.65\columnwidth}{!}{%

				\begin{tabular}{|p{1.5cm}|c|c|c|c|c|c|}

					\hline \hline
					& 1 cons & 2 cons & 3 cons & 4 cons & 5 cons   \\ \hline
					DE Hard & 4\% & 4\%  &  3\%   &  2 \%  & 2\%  \\ \hline
					ES Hard & 14\% &  10\% &    8 \%  &   7\% &  5\% \\ \hline
					PSO Hard & 15\% & 10\%  &   9\%  &   8\% &  7\%  \\ \hline

				\end{tabular}}
				\setlength{\tabcolsep}{5em}

				\label{table:feasibilityratioquad}
			\end{table}

\section{Conclusion}
In this paper, we carried out an algorithm performance comparison on each others constrained problem instances. We then analysed the features and characteristics of constraints that make them hard to solve for certain algorithm but easy for the others. It is observed that some constraint features are more contributing to problem difficulty for certain algorithms. In linear constraints, some features such as coefficient relationship, angle, local-optimum feasibility ratio and shortest distance play an important role in problem difficulty to DE and ES algorithms. Considering quadratic instances, angle does not show any relationship to problem difficulty.

By analysing how well one algorithm performs in conditions where other algorithms fail, we can derive its strengths and weaknesses over constrained problems. These results can help us to improve the efficiency of algorithm prediction model.

\section*{Acknowledgements}
Frank Neumann has been supported by ARC grants DP130104395 and DP140103400.

\begin{figure*}
	\hbox{\hspace{-18ex}
		\includegraphics[height=0.85\textwidth]{std_lin_quadArchive}}
	\caption{Box plot for standard deviation of coefficients in linear constraints with objective functions: Sphere (A), Ackley (C) and Rosenbrock (E) and quadratic constraints Sphere (B), Ackley (D) and Rosenbrock (F). Each sub figure includes hard instances (H) with 1 to 5 constraints using algorithms (a/b/c denotes a: number of constraints, b: hard instances and c: hard instances for DE/ES/PSO algorithm).}
	\label{fig:box_plot_std}

\end{figure*}

\begin{figure*}
	\hbox{\hspace{-18ex}

		\includegraphics[height=0.9\textwidth]{dis_lin_quadArchive}}
	\caption{Box plot for shortest distance feature in linear constraints with objective functions: Sphere (A), Ackley (C) and Rosenbrock (E) and quadratic constraints Sphere (B), Ackley (D) and Rosenbrock (F). Each sub figure includes hard instances (H) with 1 to 5 constraints using algorithms (a/b/c denotes a: number of constraints, b: hard instances and c: hard instances for ES/PSO/DE algorithm).}
	\label{fig:box_plot_dis}
\end{figure*}