∎

11institutetext:  1 Department of Artificial Intelligence, Key Laboratory of
Digital Protection and Intelligent Processing of Intangible Cultural Heritage
of Fujian and Taiwan, Ministry of Culture and Tourism, School of Informatics,
Xiamen University, Fujian, China, 361005.
2 11email<EMAIL_ADDRESS>
3 11email<EMAIL_ADDRESS>
4 Department of Computing, Hong Kong Polytechnic University, Hong Kong.
11email<EMAIL_ADDRESS>

# Solving Large-Scale Multi-Objective Optimization via Probabilistic
Prediction Model ††thanks: Min Jiang is the corresponding author.
This work was supported by the National Natural Science Foundation of China
(No.61673328).

Haokai Hong 1 Kai Ye 1
Min Jiang 1,2${}^{(\textrm{{\char 0\relax}})}$ Donglin Cao
1,3${}^{(\textrm{{\char 0\relax}})}$ Kay Chen Tan 4

(Received: date / Accepted: date)

###### Abstract

The main feature of large-scale multi-objective optimization problems (LSMOP)
is to optimize multiple conflicting objectives while considering thousands of
decision variables at the same time. An efficient LSMOP algorithm should have
the ability to escape the local optimal solution from the huge search space
and find the global optimal. Most of the current researches focus on how to
deal with decision variables. However, due to the large number of decision
variables, it is easy to lead to high computational cost. Maintaining the
diversity of the population is one of the effective ways to improve search
efficiency. In this paper, we propose a probabilistic prediction model based
on trend prediction model and generating-filtering strategy, called LT-PPM, to
tackle the LSMOP. The proposed method enhances the diversity of the population
through importance sampling. At the same time, due to the adoption of an
individual-based evolution mechanism, the computational cost of the proposed
method is independent of the number of decision variables, thus avoiding the
problem of exponential growth of the search space. We compared the proposed
algorithm with several state-of-the-art algorithms for different benchmark
functions. The experimental results and complexity analysis have demonstrated
that the proposed algorithm has significant improvement in terms of its
performance and computational efficiency in large-scale multi-objective
optimization.

###### Keywords:

Evolutionary multi-objective optimization Large-scale optimization
Probabilistic prediction model Trend prediction model

##  1 Introduction

Many optimization problems in the real world involve hundreds or even
thousands of decision variables and multiple conflicting goals. These problems
are defined as large-scale multi-objective optimization problems (LSMOP) RN104
; RN103 . Large-scale multi-objective optimization problems have emerged from
many real-world applications RN231 . For example, in the problem of designing
protein molecular structure, thousands of decision variables are considered.
Besides, in a routing system with large number of network nodes, the
distribution of transmission capacity of each node is regarded as a decision
variable, whereas the energy consumption, robustness and stability of the
network system constitute a large-scale multi-objectives optimization problem
7533424 . Since the size of the search space has an exponential relationship
with the number of decision variables, which leads to the problem of
dimensionality explosion in LSMOP RN105 , finding the optimal solution in
LSMOP is more challenging than ordinary multi-objective optimization
problems(MOP) 5415586 . Therefore, an excellent large-scale multi-objective
optimization algorithm (LSMOA) should overcome this problem through searching
in the decision space effectively, which is essential for solving LSMOP.

The fact that LSMOP is more difficult to solve than MOP with a small number of
decision variables makes it challenging for the multi-objective evolutionary
algorithm (MOEA) to effectively explore the search space. Generally speaking,
MOEA is tend to converge to a local optimum prematurely, or it may converge to
a relatively large areaRN121 . Therefore, it has been proved through
experience that most existing MOEAs cannot solve LSMOPs RN122 ; RN82 . In
fact, most existing MOEAs improve the effectiveness of solving MOP by
enhancing environmental selection (that is, the strategy of selecting
solutions for the next generation from the current population) RN122 ; RN124 ;
RN125 . However, the search efficiency is not sufficient for finding the
optimal solution in a large search space. Thus, it is important to find a
LSMOA that can search for the optimal solution in the space of LSMOPs.

In order to solve LSMOP, various LSMOAs have been presented in recent years
2008Large ; 8720021 ; 9047876 ; 8765790 . These proposed LSMOAs can be roughly
divided into the following three categories. The first category of LSMOAs is
benefits from the decision variable analysis. A large-scale evolutionary
algorithm based on the decision variable clustering (LMEA) RN82  is
representative in this category, which divides the decision variables into two
categories via clustering method. Then, by focusing on convergence and
diversity respectively, two different strategies are used to optimize the
variables related to convergence and diversity. The second category of LSMOAs
is based on the decision variable grouping. This type of LSMOA divides the
decision variables into several groups, and then optimizes each group of
decision variables. For example, C-CGDE3 RN84  maintains several subsets of
equal-length decision variables obtained by variable grouping as
subpopulations. The third category is based on problem conversion. A
representative framework of this category is called Large-scale Multi-
Objective Optimization Framework (LSMOF) RN88 . In LSMOF, the original LSMOP
is converted into a low-dimensional single-objective optimization problem with
some direction vectors and weight variables, the method aimed to guide the
population towards the optimal solution.

In the methods of decision variable analysis and decision variable grouping,
due to that fact that the algorithm relies on the research on decision
variables, the algorithm complexity is relatively high. However solving LSMOP
should not only consider how to escape from local optimum, but also how to
reduce the impact of the number of decision variables on the complexity of the
algorithm 10.5555/3298239.3298367 . Inspired by the advantages of Trend
Prediction Model (TPM), we found that the Trend Prediction Model can be used
to predict the optimal solution. The TPM is used to estimate and predict the
probability of potential particles at each position in the next generation,
thereby increasing the uncertainty of particle movement and achieving better
prediction results. With the advantages of trend prediction model (TPM), the
method not only maintains the diversity of the population 9122031 , but also
keeps the computational complexity independent of decision variables. In the
specific implementation, TPM was applied to sample individuals from the parent
and expand into the offspring population.

In this paper, we integrate the Trend Prediction Model into a probabilistic
prediction model to obtain a novel LSMOA algorithm. This method is a
probabilistic prediction model for solving large-scale multi-objective
optimization based on trend prediction model (LT-PPM). Our proposed algorithm
selects individuals and expands them from the previous generation to achieve
the purpose of expanding and generating a new offspring population via TPM,
which ensures the diversity of the offspring population. The contributions of
this work are summarized as follows.

  1. 1.

The proposed algorithm framework is an effective tool to solve the problem of
large decision variable space and overcome the tendency to fall into a local
optimum in LSMOP. Our method uses the trend prediction model to maintain the
quality and diversity of the population simultaneously to escape from the
local optimum, which provides a novel idea to solve LSMOP.

  2. 2.

Taking into account the high dimensionality of the decision variables of the
LSMOP problem, the complexity of our proposed method is not directly related
to the scale of the decision variables, but only related to the number of
objectives and populations. While escaping from the local optimum, the
influence of large-scale decision variables on the computational complexity is
also eliminated in our method.

##  2 Preliminary Studies and Related Works

In this section, the definition of large-scale multi-objective optimization
problems is presented. And then existing methods to solve LSMOPs are
introduced.

###  2.1 Large-scale Multi-objective Optimization

The mathematical form of LSMOPs is as follows:

$\displaystyle\mathbf{Minimize}\ F(x)=$ $\displaystyle\left\langle
f_{1}(x),f_{2}(x),...,f_{m}(x)\right\rangle$ (1) $\displaystyle s.t.\
x\in\Omega$

where $x=(x_{1},x_{2},\dots,x_{n})$ is the $n$-dimensional decision vector,
$F=(f_{1}(x),f_{2}(x),...,f_{m}(x))$ is the $m$-dimensional objective vector.

It should be noticed that the dimension of decision vector $n$ is more than
100 in LSMOP RN85 . When dealing with multi-objective optimization problems,
evolutionary algorithms (EAs) are usually based on individual evolution
strategies, genetic operators are used to evolve individuals, so that EAs have
good global convergence and versatility. The evolutionary framework we use in
this paper is an emerging branch in this field called the Probabilistic
Prediction Model (PPM) JIANG2018203 ; 8285444 , PPM promots the evolution of
the entire population by establishing a probability model on the Pareto-
optimal Front.

###  2.2 Related Works

In this part, we will introduce some state-of-art large-scale multi-objective
optimization algorithms. Much progress has been made in solving LSMOPs. Most
of existing methods can be divided into three categories: based on decision
variable analysis, based on decision variable grouping and based on problem
transformation.

The first category is based on the analysis of decision variables. Zhang et
al. RN82  proposed a large-scale evolutionary algorithm (LMEA) based on
clustering of decision variables. Another representative of this type of
algorithm is MOEA/DVA RN106 , in which the original LSMOP is decomposed into
many simpler sub-MOPs. Then, the decision variables in each sub-problem are
optimized as an independent sub-component. The disadvantage of decision
variable analysis-based LSMOAs is that the connection between decision
variables could be incorrectly determined to update the solution, which may
lead to local optimization RN82 . To solve large-scale MOPs that pose a great
challenge to the balance between diversification and convergence by
considering a correlation between position and distance functions, a novel
indicator-based algorithm with an enhanced diversification mechanism was
developed by Hong et al. 8533425 .

The second type of LSMOAs is decision variable grouping-based framework RN83 .
The method proposed by Antonio et al. RN84  maintained several independent
subgroups. Each subgroup is a subset of the same-length decision variables
obtained by variable grouping (for example, random grouping, linear grouping,
ordered grouping or differential grouping). All sub-populations work together
to optimize LSMOP in a divide-and-conquer manner. Zille et al. RN89  proposed
a weighted optimization framework (WOF) for solving LSMOP. Decision variables
are divided into many groups, and each group is assigned a weight variable.
Then, in the same group, the weight variable is regarded as a sub-problem of a
sub-space from the original decision space. Liu et al. LIU2020100684  proposed
an algorithm in which the random dynamic grouping was used to substitute
ordered grouping to improve the performance of the WOF framework.

In this category, there is a novel idea that based on competitive optimizer.
Tian et al. RN85  proposed LMOCSO based on a competitive population optimizer
to solve LSMOP. The method used a two-stage strategy to update the position of
the particles. Another idea of this category is based on collaboration. Bergh
et al. RN87  proposed a method to transform particle swarm optimization
algorithm into a collaborative framework. It divides the search space into
low-dimensional subspaces to overcome the exponential increase in difficulty
and ensures that it can search every possible area of the search space. And
Shang et al. RN86  proposed the idea of cooperative co-evolution, which was
used to deal with large-scale multi-target capability routing problems.
Another cooperative coevolutionary algorithm was proposed by Li et al.
10.1145/3205651.3208250 , a fast interdependency identification grouping
method was utilized for large-scale multi-objective problems.

However, for LSMOA based on the grouping of decision variables, two related
variables can be divided into different groups, which may lead to local
optimization RN83 .

The third category is based on problem transformation. In order to enhance the
pressure of convergence, one intuitive solution is to directly convert the
definition of traditional Pareto dominance, such as $\theta$-dominance RN116 ,
$L$-optimality RN113 , $\epsilon$-dominance RN112 , preference order ranking
RN115  and fuzzy dominance RN114 . To enhance the convergence performance of
evolutionary algorithm in solving large-scale multi-objective problems, Qin et
al. 9367299  proposed a method for generating guiding solutions along search
directions defined by a few chosen individuals in the current parent
population. Another kind of idea that belongs to this category is to combine
traditional dominance with other convergence-related metrics. MOEAs that fall
into this category include Grid-based Evolutionary Algorithm (GrEA) RN118 ,
NSGA-II based on substitute distance assignment RN117 , many-Objective
Evolutionary Algorithm based on Directional Diversity and Favorable
Convergence (MaOEA-DDFC) RN120 , Preference-Inspired Co-Evolutionary Algorithm
(PICEA-g) RN119 , and Knee point driven Evolutionary Algorithm (KnEA) RN109 .

In addition to the above-mentioned convergence-related problem transformation,
for LSMOP, He et al. RN88  introduced a general framework called the large-
scale multi-objective optimization framework (LSMOF). LSMOF reconstructs the
original LSMOP into a low-dimensional single-objective optimization problem
with some direction vectors and weight variables, aiming to lead the
population to POS. Recently, He et al. RN90  designed an evolutionary multi-
objective optimization (GMOEA) driven by GANs, which produces offspring from
GANs.

As the number of decision variables increases linearly, the complexity (as
well as volume) of the search space will increase exponentially, and the
methods of analysis or grouping of decision variables presented above cannot
tackle the curse of dimensionality. Therefore, in this work, we focus on how
to use TPM to maintain population diversity to avoid falling into local
optimum while avoiding the curse of dimensionality caused by large-scale
decision variables to better solve LSMOP.

##  3 A Probabilistic Prediction Model Based on Trend Prediction Model

###  3.1 Overview

Based on the Generating-Filtering strategy, LT-PPM combines the importance
sampling method with the non-parametric density estimation function to propose
a different selection method for candidate solution. Specifically, LT-PPM uses
a kernel density estimation method, a non-parametric probability density
estimation method to assign a sparseness to each individual, and then uses
importance sampling to select subsequent solutions. This method of selecting
candidate solutions aims to ensure the balance between exploitation and
exploration in terms of time and space: on the one hand, the sparseness will
change with generation, on the other hand, sparseness itself will also guide
the method to achieve the balance in the large-scale search space.

We first give the definition of particles, and then in the remainder of this
chapter we will introduce importance sampling based on kernel density
estimation, trend prediction model, LT-PPM and the complexity analysis of the
algorithm.

######  Definition 1 (Particle)

For an optimization problem with $n$ optimization goals, the particle is a
two-tuple $p=(x,\boldsymbol{v})$, where $x$ represents a solution of the
multi-objective optimization problem, and $\boldsymbol{v}$ is a unit vector in
$n$-dimensional space representing the direction of the particle.

Note that the direction vector of the first generation population is randomly
initialized. In the subsequent offspring population, the direction vector of
each individual is the direction of movement when their parent particle
produce them.

###  3.2 Importance Sampling Based on Kernel Density Estimation

Importance Sampling is defined as follows,

######  Definition 2 (Importance Sampling)

When calculating the expectation of the statistic $f(x)$ on the original
probability distribution $p(x)$, one can calculate the expectation of
$f(x)\cdot p(x)/q(s)$ on a certain importance probability distribution $q(x)$,

$\displaystyle\mathbf{E}_{p}[f(x)]$ (2) $\displaystyle=$
$\displaystyle\int_{x}p(x)\cdot f(x)dx$ $\displaystyle=$
$\displaystyle\int_{x}q(x)\cdot f(x)\frac{p(x)}{q(x)}dx$ $\displaystyle=$
$\displaystyle\mathbf{E}_{q}[f(x)\frac{p(x)}{q(x)}]$

Assuming that the current population is
$P=\\{x_{1},\dots,x_{n}\\}$, and corresponding particles are
$P^{\prime}=\\{p_{1},\dots,p_{n}\\}$. The probability density of the
population in different regions reflects the degree of density in the search
space: where the probability density is high, the population is denser.
Therefore, the density of particle $p_{k}$ can be defined as its kernel
density estimate,

$\displaystyle\mathbf{Density}(p_{k})=\frac{1}{nh}\sum_{i=1}^{n}\kappa(\frac{p_{k}-p_{i}}{h})$
(3)

where $\kappa(\cdot)$ is the kernel function, and the parameter $h$ is used to
adjust the width of the kernel density curve: when $h$ is larger, the
distribution curve is smoother, when $h$ is smaller, the distribution curve is
steeper.

The Sparseness of particle $p_{k}$ is defined as the reciprocal of Density,

$\displaystyle\mathbf{Sparseness}(p_{k})=\frac{1}{\mathbf{Density}(p_{k})}$
(4)

The probability density of the population in different regions reflects the
degree of density in the search space: where the probability density is high,
the population is denser. Therefore, the density of particle $p_{k}$ can be
defined as its kernel density estimate, The Sparseness of particle $p_{k}$ is
defined as the reciprocal of Density. Based on the Sparseness of each
particle, Importance Sampling for our method can be defined. Since each
particle is randomly selected, particle $p_{k}$ obeys a uniform distribution
$U(p_{k})=1/n$. Then we have importance distribution,

$\displaystyle
IS(p_{k})=\frac{\mathbf{sparseness}(p_{k})}{\sum_{i}\mathbf{sparseness}(p_{i})}$
(5)

The number of offspring that the particle $p_{k}$ will produce is proportional
to $S(p_{k})$. The expectation value of sparseness describes the sampling
process, and its calculation formula is: $E_{U}[S(p_{k})]$. However, in MOPs,
it is not easy to determine how many offspring the selected particle should
produce, so we do not want to directly use this method for sampling. The
original sampling process can be rewritten as,

$\displaystyle\mathbf{E}_{U}[S(p_{k})]=\mathbf{E}_{S}[S(p_{k})\cdot\omega(p_{k})]=\mathbf{E}_{S}[U(p_{k})]$
(6)

where $\omega(p_{k})=\frac{U(p_{k})}{S(p_{k})}$.

This means that we can select particles by the importance distribution $S$,
and then produce the same number of offspring. Algorithm 1 summarizes the
above sampling process. In each iteration, first, the sparseness of each
particle is calculated, and then the importance distribution $S$ is
constructed, then an individual is sampled on the distribution $S$, and the
individual and its velocity are returned. This sampling method does not rely
on the analysis of decision variables, avoiding the influence of large-scale
decision variables on computational complexity.

Input: $POS$ (Pareto-optimal Set), $h$ (bandwidth)

Output: $p_{i}$ (sampled particle), $x_{k}$(solution of $p_{i}$),
$\boldsymbol{v_{k}}$(velocity of $p_{i}$)

1 $n=$ size of $POS$ ;

2 for  _$p_{i}$ in $POS$_ do

3
$\mathbf{Density}(p_{k})=\frac{1}{nh}\sum_{i=1}^{n}\kappa(\frac{p_{k}-p_{i}}{h})$
;

4 $\mathbf{Sparseness}(p_{k})=\frac{1}{\mathbf{Density}(p_{k})}$ ;

5

6for  _$x_{i}$ in $POS$_ do

7 Constructe $IS(x_{i})$ according to Equation 5 ;

8

9Sampled $x_{k}$ based on $IS(x_{i})$ ;

return $x_{k}$, $\boldsymbol{v_{k}}$

Algorithm 1 ISample($POS$)

###  3.3 Trend Prediction Model

After selecting the particles to be expanded via above selection process, it
is necessary to sample on a specific distribution represented by selected
particles to generate new individuals. Physically, if an object is not
affected by external forces, it will maintain its original state of motion.
Therefore, if a particle moves in the direction $\boldsymbol{v}$, the
displacement direction $\boldsymbol{s}$ of the particle at the next moment is
equal to $\boldsymbol{v}$.

In order to effectively move to the POF, the randomness is added to the
movement of the particle, so that the direction of movement only provides
information on the trend of the movement instead of directly determining the
place at the next moment. Such trend information can well guide the particles
to update, thereby approaching the POF faster and more accurately.

######  Definition 3 (Trend Prediction Model, TPM)

For any given direction vector $\boldsymbol{v}\in\mathbb{R}^{n}$ and distance
measurement function $d(\cdot,\cdot)$. A probability distribution $P$ on
$D\subset\mathbb{R}^{n}$ is called as a trend prediction model (TPM).

For the two vectors $\boldsymbol{t_{1}}$ and $\boldsymbol{t_{2}}$ in the
model, if it satisfies:

$\displaystyle\forall\boldsymbol{t_{1}},\boldsymbol{t_{2}}\in
D,d(\boldsymbol{v},\boldsymbol{t_{1}})\leq
d(\boldsymbol{v},\boldsymbol{t_{2}})\to P(\boldsymbol{t_{1}})\geq
P(\boldsymbol{t_{2}})$ (7)

then it means that the vector $t_{1}$ is closer to the direction vector $v$
than $t_{2}$, and $t_{1}$ should have a higher probability of occurrence.

Input: $\boldsymbol{v}\in\mathbb{R}^{n}$ (Unit direction vector)

Output: $\boldsymbol{u}$ (Random direction vector)

1 $n=$ Dimension of $\boldsymbol{v}$;

2 $t=$ The subscript of the first non-zero vector in $\boldsymbol{v}$;

3 for  _$i=1,2,\dots,n-1$_ do

4 if  _$i=t$_ then

5 $\boldsymbol{\beta_{i}}$ = The vector that only the $n$th element is $1$ and
other elements are $0$;

6 else

7 $\boldsymbol{\beta_{i}}$ = The vector that only the $i$th element is $1$ and
other elements are $0$;

8

9$R$ = Apply $\mathbf{Schmidt\ orthogonalization}$ to the matrix
$[\boldsymbol{v},\boldsymbol{\beta_{1}},\dots,\boldsymbol{\beta_{n-1}}]$;

10 $\boldsymbol{u}$ = Sampled from the $n-1$ dimensional unit hyperspherical
coordinate system ;

11 $\boldsymbol{u}=R^{-1}\cdot[\cos(\theta)\boldsymbol{u}]^{T}$ ;

return $\boldsymbol{u}$

Algorithm 2 randpick($\boldsymbol{v},\theta$)

In the trend prediction model, particle movement has a moving trend, which is
indicated by the unit vector $\boldsymbol{v}$. The method aims to design the
TPM distribution so that the closer the sampling vector to $\boldsymbol{v}$,
the higher the sampling probability. For this, $\mathbf{TPM}_{h}$ is
introduced. Given the distance measurement method $d$ is the $\cos$ distance,
the domain $D$ is the unit hypersphere, and sampled from the Gaussian
distribution with the expectation of 0 and the variance of $1/h$, which is the
bandwidth. Since the domain is a unit hypersphere, $\mathbf{TPM}_{h}$ is
actually a probability distribution on a direction vector. In this
distribution, directions close to the vector $\boldsymbol{v}$ have a higher
probability of being selected (that is, the angle with $\boldsymbol{v}$ is
small), otherwise the probability is low.

After sampling an angle $\theta$ on the $\mathbf{TPM}_{h}$ system, it is
necessary to randomly select one among all the vectors whose angle is $\theta$
with $\boldsymbol{v}$. In order to sample all vectors whose angle is $\theta$
with $\boldsymbol{v}$, we can construct a rotation matrix $R$ via Schmidt
Orthogonalization to rotate the vector $\boldsymbol{v}$ to the $X$ axis. Then,
randomly select a vector on the $n-1$ dimensional hypersphere. Finally, the
sampled vector is rotated back to the original coordinate space through the
inverse transformation $R^{(-1)}$.

Algorithm 2 is the process of generating offspring direction vectors according
to the parent individual by importance sampling. The complete trend prediction
model is given in Algorithm 3, where ISample (POS) is Algorithm 1. In the
specific strategy of generating offspring, after sampling the direction vector
by Algorithm 2, it is necessary to continue to determine the particle’s
movement step length. The actual step length is sampled from a normal
distribution with a expectation of 0 and a variance of $h$.

Input: $POS$ (Pareto-Optimal Set), $h$ (bandwidth), $N$ (number of samples)

Output: $S$(Sampled set)

1 $S=\emptyset$ ;

2 for  _$i=1,2,\dots,N$_ do

3 $<p,v>=\mathbf{ISample}(POS)$ ;

4 $\theta$ = Sampled from $\mathbf{Norm}(0,1/h)$ ;

5 $u=\mathbf{randpick}(v,\theta)$ ;

6 $\lambda=$ Sampled from $\mathbf{Norm}(0,h)$ ;

7 $\lambda=\lceil\lambda\rceil$ ;

8 $S=S\cup\\{p+\lambda u\\}$ ;

9

10return $S$ ;

Algorithm 3 TPM($POS,h,N$) (Trend Prediction Model)

###  3.4 Large scale TPM based Probabilistic Prediction Model (LT-PPM)

Combining the importance sampling and trend prediction model, the complete
algorithm implementation of LT-PPM is given in Algorithm 4. The key steps of
the algorithm are:

Generating: Select the particles to be expanded through importance sampling
based on kernel density estimation, and then predict the next generation
population through the trend prediction model;

Filtering: Screen out outstanding individuals from the parent population and
offspring population.

Input: $N$(size of restricted population), $P$(size of population), $e$ (the
maximum evaluations), $r$(attenuation coefficient)

Output: $POS(the\ \mathbf{POS})$

1 Random initial $P$ population $POS$ ;

2 while  _Used evaluations $\leq e$_ do

3 $S=sample(POS,1/h,P)$ ;

4 Update $POS$ and based on $S$ ;

5 Remove the densest individuals repeatedly until $|POS|\leq N$ ;

6 $h=h\cdot r$;

7

8return $POS$ ;

Algorithm 4 Large scale TPM based PPM

The algorithm framework of LT-PPM is briefly illustrated in conjunction with
Figures 1(a) to 1(d). In a certain generation, the LT-PPM algorithm maintains
a population and calculates the corresponding POS and POF (the solid line in
Figure 1(a) represents the POF of this generation, and the dashed line
represents the real POF). Then use the non-parametric kernel density
estimation method to estimate the probability density of the distribution
corresponding to the POS, and use the importance sampling method to select a
particle. The principle is that particles with higher sparseness will have a
higher probability of being selected (Figure 1(b)), while particles in dense
particles have a lower probability of being selected. After that, apply the
trend prediction model to the selected particles to obtain a series of new
solutions (Figure 1(c), where gray dots represent the newly generated
population). Finally, update POS and POF according to the newly generated
individuals (Figure 1(d)).

Figure 1: Algorithm Framework of LT-PPM

###  3.5 Complexity Analysis

Assuming a $n$ objectives LSMOP with $m$ population, then the time to
calculate all objectives is $\mathcal{O}(n)$. In our algorithm, $m$
individuals are sampled in each iteration, and each individual requires
$\mathcal{O}(m)$ time to sample, therefore the time complexity of the sampling
process in each iteration is $\mathcal{O}(m^{2})$. Suppose the evaluation is
$e$, the actual interation times is $e/m$, so the total complexity of the
algorithm is $\mathcal{O}(nme)$.

The complexity of the algorithm is proportional to the square of the number of
population and the number of objectives, and has no direct relation with the
scale of decision variables of the problem. Therefore, our proposed LT-PPM can
avoid the influence of large-scale decision variables on the computational
complexity in solving large-scale multi-objective optimization problems.

##  4 Experiments and Analysis

###  4.1 Algorithms in Comparison and Test Problems

The proposed LT-PPM was experimented on large-scale multi-objective problems
LSMOP1 - LSMOP9 RN96  and compared with several state-of-the-art algorithms,
including LMEA RN82 , LMOCSO RN85 , NSGAIII RN110  and RMMEDA RN111 .
According to RN92 , the parameters of LMEA, RMMEDA and LMOCSO algorithms are
set to default values. For all test problems, the number of optimization
objectives is set to 3. The remaining parameters are set as follows,

  1. 1.

Population size: The population size $N$ is set to 300.

  2. 2.

Termination condition: The maximum evaluation times $e$ of all compared LSMOAs
is set to 100000.

  3. 3.

Parameter of LT-PPM: The bandwidth attenuation coefficient $r$ in algorithm 4
is set to 0.9.

###  4.2 Performance Indicators

  1. 1.

Schott’s Spacing (SP) RN97 : It measures the uniformity of the solution found
by the algorithm. The smaller the SP value, the more uniform the distribution
of solutions obtained by the algorithm. It is calculated by the following
formula:

$\displaystyle
SP=\sqrt{\frac{1}{n_{\text{POF}^{*}}-1}\cdot\sum_{i=1}^{n_{\text{POF}^{*}}}(E_{i}-\bar{E})^{2}}$
(8)

where $\text{POF}^{*}$ is the true POF of the multi-objective optimization
problem, $E_{i}$ represents the Euclidean distance between the $i$-th solution
in $\text{POF}^{*}$ and its closest solution, and $\bar{E}$ represents the
mean value of $E_{i}$.

  2. 2.

Inverted Generational Distance (IGD) RN98 : IGD is a metric used to quantify
the convergence of the solution obtained by the multi-objective optimization
algorithm. When the IGD value is small, the convergence of the solution is
improved. IGD is defined as

$\displaystyle
IGD(\text{POF}^{*},\text{POF})=\frac{1}{n}\cdot\sum_{p^{*}\in\text{POF}^{*}}\min_{p\in\text{POF}}||p^{*}-p||^{2}$
(9)

where $\text{POF}^{*}$ is the true POF of the multi-objective optimization
problem, POF is the approximate set of POF obtained by the multi-objective
optimization algorithm, and $n$ is the number of solutions in
$\text{POF}^{*}$.

###  4.3 Performance on LSMOP Problems and Discussion

Table 1: Mean Values of SP Metric Obtained by Compared Algorithms Over LSMOP Problems With 3-Objectives Problem | Des. | NSGAIII | RMMEDA | LMEA | LMOCSO | LT-PPM
---|---|---|---|---|---|---
| 1000 | 4.37E-01 | 4.36E-01 | 4.72E-01 | 4.75E-01 | 4.15E-01
LSMOP1 | 2000 | 4.43E-01 | 4.28E-01 | 4.33E-01 | 5.16E-01 | 4.93E-01
| 5000 | 4.35E-01 | 4.32E-01 | 4.19E-01 | 5.33E-01 | 3.15E-01
| 1000 | 1.68E-02 | 4.80E-01 | 5.34E-01 | 6.19E-02 | 3.10E-01
LSMOP2 | 2000 | 1.37E-02 | 4.81E-01 | 5.22E-01 | 5.27E-02 | 3.10E-01
| 5000 | 1.47E-02 | 5.44E-01 | 5.35E-01 | 5.24E-02 | 3.16E-01
| 1000 | 8.62E-01 | 9.27E-01 | 1.07E+00 | 6.67E-01 | 1.76E-01
LSMOP3 | 2000 | 8.77E-01 | 9.72E-01 | 1.02E+00 | 8.36E-01 | 2.09E-01
| 5000 | 9.23E-01 | 9.03E-01 | 1.09E+00 | 9.86E-01 | 3.66E-01
| 1000 | 5.24E-02 | 4.81E-01 | 5.12E-01 | 1.11E-01 | 2.97E-01
LSMOP4 | 2000 | 2.63E-02 | 5.13E-01 | 5.22E-01 | 7.39E-02 | 3.10E-01
| 5000 | 1.63E-02 | 5.14E-01 | 5.19E-01 | 5.62E-02 | 2.99E-01
| 1000 | 6.61E-01 | 7.61E-01 | 7.25E-01 | 7.13E-01 | 4.39E-01
LSMOP5 | 2000 | 5.75E-01 | 7.04E-01 | 7.74E-01 | 6.58E-01 | 4.10E-01
| 5000 | 5.53E-01 | 6.86E-01 | 7.68E-01 | 7.04E-01 | 3.70E-01
| 1000 | 9.77E-01 | 1.02E+00 | 1.15E+00 | 1.30E+00 | 3.08E-01
LSMOP6 | 2000 | 9.44E-01 | 1.03E+00 | 1.15E+00 | 1.24E+00 | 2.58E-01
| 5000 | 9.57E-01 | 1.15E+00 | 1.17E+00 | 9.31E-01 | 3.07E-01
| 1000 | 8.59E-01 | 8.73E-01 | 9.46E-01 | 3.57E-01 | 2.86E-01
LSMOP7 | 2000 | 5.97E-01 | 9.11E-01 | 9.11E-01 | 1.02E+00 | 3.13E-01
| 5000 | 5.21E-01 | 9.37E-01 | 8.99E-01 | 3.91E-01 | 3.44E-01
| 1000 | 5.24E-01 | 7.39E-01 | 7.23E-01 | 9.56E-01 | 3.39E-01
LSMOP8 | 2000 | 4.54E-01 | 7.53E-01 | 7.63E-01 | 8.58E-01 | 3.29E-01
| 5000 | 3.89E-01 | 7.07E-01 | 7.42E-01 | 8.99E-01 | 3.19E-01
| 1000 | 8.58E-01 | 7.54E-01 | 9.46E-01 | 8.93E-01 | 3.66E-01
LSMOP9 | 2000 | 8.63E-01 | 8.16E-01 | 9.05E-01 | 9.67E-01 | 3.87E-01
| 5000 | 8.50E-01 | 8.26E-01 | 9.17E-01 | 9.77E-01 | 5.32E-01
Table 2: Mean Values of IGD Metric Obtained by Compared Algorithms Over LSMOP Problems With 3-Objectives Problem | Des. | NSGAIII | RMMEDA | LMEA | LMOCSO | LT-PPM
---|---|---|---|---|---|---
| 1000 | 5.4E+00 | 8.8E+00 | 1.1E+01 | 1.3E+00 | 6.5E-01
LSMOP1 | 2000 | 7.3E+00 | 7.8E+00 | 1.1E+01 | 1.5E+00 | 5.9E-01
| 5000 | 9.8E+00 | 8.6E+00 | 1.1E+01 | 1.4E+00 | 7.5E-01
| 1000 | 3.1E-02 | 3.4E-02 | 3.4E-02 | 2.6E-02 | 3.9E-02
LSMOP2 | 2000 | 2.0E-02 | 2.5E-02 | 2.4E-02 | 1.9E-02 | 3.1E-02
| 5000 | 1.5E-02 | 2.0E-02 | 2.0E-02 | 1.4E-02 | 2.7E-02
| 1000 | 1.4E+01 | 1.6E+01 | 5.5E+01 | 9.9E+00 | 4.9E+00
LSMOP3 | 2000 | 1.7E+01 | 1.6E+01 | 4.2E+01 | 1.1E+01 | 4.7E+00
| 5000 | 1.9E+01 | 1.7E+01 | 2.1E+02 | 1.1E+01 | 5.1E+00
| 1000 | 1.1E-01 | 1.1E-01 | 1.1E-01 | 8.3E-02 | 1.1E-01
LSMOP4 | 2000 | 5.8E-02 | 6.3E-02 | 6.2E-02 | 4.8E-02 | 6.4E-02
| 5000 | 2.8E-02 | 3.2E-02 | 3.2E-02 | 2.5E-02 | 3.7E-02
| 1000 | 9.5E+00 | 1.1E+01 | 1.9E+01 | 3.0E+00 | 6.4E-01
LSMOP5 | 2000 | 1.2E+01 | 1.2E+01 | 1.9E+01 | 3.4E+00 | 4.8E-01
| 5000 | 1.8E+01 | 1.3E+01 | 2.0E+01 | 3.3E+00 | 7.6E-01
| 1000 | 3.7E+03 | 7.2E+03 | 3.1E+04 | 3.5E+02 | 1.9E+00
LSMOP6 | 2000 | 8.8E+03 | 1.2E+04 | 3.5E+04 | 5.0E+02 | 1.0E+00
| 5000 | 1.8E+04 | 5.8E+03 | 3.2E+04 | 5.6E+02 | 8.7E+00
| 1000 | 1.1E+00 | 1.1E+00 | 5.4E+02 | 8.2E-01 | 2.8E+00
LSMOP7 | 2000 | 1.0E+00 | 1.0E+00 | 1.4E+02 | 7.2E-01 | 4.0E+00
| 5000 | 9.7E-01 | 9.7E-01 | 1.2E+02 | 6.5E-01 | 2.1E+01
| 1000 | 5.8E-01 | 5.8E-01 | 5.9E-01 | 5.0E-01 | 2.0E-01
LSMOP8 | 2000 | 5.6E-01 | 5.8E-01 | 5.6E-01 | 4.9E-01 | 1.9E-01
| 5000 | 5.6E-01 | 5.5E-01 | 5.6E-01 | 4.9E-01 | 1.9E-01
| 1000 | 2.2E+01 | 4.1E+01 | 1.4E+02 | 2.4E+00 | 1.4E+00
LSMOP9 | 2000 | 4.3E+01 | 7.3E+01 | 1.4E+02 | 8.8E+00 | 1.4E+00
| 5000 | 6.9E+01 | 8.7E+01 | 1.5E+02 | 4.2E+01 | 1.4E+00

Table 1 and Table 2 show the statistical results of the average SP and IGD
values for 10 runs respectively, where the best result on each test instance
is shown in bold.

From Table 1 and Table 2, it obviously that the performance of the proposed
LT-PPM is better than the other five comparative convergence algorithms. Table
1 lists the SP values of the six comparison algorithms. LT-PPM performed best
in 20 out of 27 test instances, followed by NSGAIII with 6 best results, and
RMMEDA with 1 best result. Specifically, under all decision variable settings,
LT-PPM performed better on LSMOP3, LSMOP5, LSMOP6, LSMOP7, LSMOP8 and LSMOP9.
In other test instances, LT-PPM is slightly lower than the corresponding best
performance algorithm.

In Table 2, LT-PPM got 18 of the 27 best results, LMOCSO got 9 best results,
and the other algorithms did not get the best results. Specifically, under all
decision variable settings, LT-PPM performs better on LSMOP1, LSMOP3, LSMOP5,
LSMOP6, LSMOP8 and LSMOP9. LMOCSO achieved better convergence on problems
LSMOP2, LSMOP4 and LSMOP7.

The above experimental results show that for most test instances, LT-PPM can
obtain a set of solutions with good convergence and diversity. It can be seen
that the proposed LT-PPM can obtain competitive performance on the LSMOP test
problem, which confirms the effectiveness of LT-PPM in solving the challenging
large-scale decision variable problem.

###  4.4 Experiment on Time Complexity and Discussion

Figure 2: Experiment on Time Complexity

Figure 2 compares the CPU running time on the LSMOP1 problem of the algorithms
in the experiment, where the abscissa is the different decision variables, and
the ordinate is the CPU time. It is worth noting that when the decision
variable is greater than 500, the running time of RM-MEDA is greatly
increased.

It can be seen from the results that the time consumption of our proposed
algorithm is not the lowest when there are fewer decision variables. However,
as the decision variables increase, the running time of the proposed algorithm
grows slower than other algorithms. When it reaches 5000 dimensions, the time
consumption of our algorithm is the lowest.

The results of the experiment prove that our algorithm is not directly related
to the decision variables. This allows the algorithm to have better
convergence and diversity when solving large-scale multi-objective
optimization problems, while also avoiding the impact of large-scale decision
variables on computational complexity.

##  5 Conclusion and Future Works

This paper proposes an algorithm, called LT-PPM, for solving LSMOP. The TPM
component of the proposed method enables us to effectively use the heuristic
information of the individual motion process to predict the POF for next
generation. At the same time, the importance sampling technique makes the
algorithm independent of the numbers of decision variables, so that the
algorithm can explore faster in the large-scale search space. Simultaneously,
the importance sampling and the TPM model are used to adjust the sparsity of
the population, thereby balancing the development and exploration of the
algorithm and escaping from the local optimum.

In LT-PPM algorithm, the general modeling-sampling process is not used.
Instead, the Generating-Filtering strategy is adopted, which enables the
algorithm to focus on using heuristic information to generate individuals and
screen out dominant populations. Our method is compared with several other
large-scale multi-objective optimization algorithms on 9 LSMOPs. The
experimental results show that LT-PPM has good performance in convergence,
distribution and robustness.

This paper presents a preliminary work on the application of trend prediction
models to large-scale evolutionary calculations. Future work may take the
following possible directions. The method has some details that can be
improved, for large-scale problems with complex search spaces, such as
constraints and discontinuous Pareto optimal regions, the proposed framework
still has a lot of room for improvement. In addition, a wealth of advanced
machine learning techniques can stimulate more innovations to solve practical
applications with large-scale decision variables.

###### Acknowledgements.

Acknowledgement This work was supported by the National Natural Science
Foundation of China (No.61673328).

## References

  * [1] A Ponsich, A. L Jaimes, and C. A. C Coello.  A survey on multiobjective evolutionary algorithms for the solution of the portfolio optimization problem and other finance and economics applications.  IEEE Transactions on Evolutionary Computation, 17(3):321–344, 2013\.
  * [2] Z. P Stanko, T Nishikawa, and S. R Paulinski.  Large-scale multi-objective optimization for the management of seawater intrusion, santa barbara, ca.  In Agu Fall Meeting, 2015.
  * [3] Wen-Jing Hong, Peng Yang, and Ke Tang.  Evolutionary computation for large-scale multi-objective optimization: A decade of progresses.  International Journal of Automation and Computing, 18(2):155–169, 2021.
  * [4] Ke Tang, Juan Wang, Xiaodong Li, and Xin Yao.  A scalable approach to capacitated arc routing problems based on hierarchical decomposition.  IEEE Transactions on Cybernetics, 47(11):3928–3940, 2017.
  * [5] Handing Wang, Licheng Jiao, Ronghua Shang, Shan He, and Fang Liu.  A memetic optimization strategy based on dimension reduction in decision space.  Evolutionary Computation, 23(1), 2015.
  * [6] Juan J. Durillo, Antonio J. Nebro, Carlos A. Coello Coello, José Garcia-Nieto, Francisco Luna, and Enrique Alba.  A study of multiobjective metaheuristics when solving parameter scalable problems.  IEEE Transactions on Evolutionary Computation, 14(4):618–635, 2010\.
  * [7] Elre T. Oldewage, Andries P. Engelbrecht, and Christopher W. Cleghorn.  The merits of velocity clamping particle swarm optimisation in high dimensional spaces.  In 2017 IEEE Symposium Series on Computational Intelligence (SSCI), pages 1–8, 2017.
  * [8] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan.  A fast and elitist multiobjective genetic algorithm: Nsga-ii.  IEEE Transactions on Evolutionary Computation, 6(2):182–197, 2002\.
  * [9] X. Zhang, Y. Tian, R. Cheng, and Y. Jin.  A decision variable clustering-based evolutionary algorithm for large-scale many-objective optimization.  IEEE Transactions on Evolutionary Computation, 22(1):97–112, 2018\.
  * [10] Y. Tian, R. Cheng, X. Zhang, F. Cheng, and Y. Jin.  An indicator-based multiobjective evolutionary algorithm with reference point adaptation for better versatility.  IEEE Transactions on Evolutionary Computation, 22(4):609–622, 2018\.
  * [11] Q. Zhang and H. Li.  Moea/d: A multiobjective evolutionary algorithm based on decomposition.  IEEE Transactions on Evolutionary Computation, 11(6):712–731, 2007\.
  * [12] Z. Yang, K. Tang, and X. Yao.  Large scale evolutionary optimization using cooperative coevolution.  Information Sciences, 178(15):2985–2999, 2008.
  * [13] Ye Tian, Xingyi Zhang, Chao Wang, and Yaochu Jin.  An evolutionary algorithm for large-scale sparse multiobjective optimization problems.  IEEE Transactions on Evolutionary Computation, 24(2):380–393, 2020\.
  * [14] Ye Tian, Chang Lu, Xingyi Zhang, Kay Chen Tan, and Yaochu Jin.  Solving large-scale multiobjective optimization problems with sparse optimal solutions via unsupervised neural networks.  IEEE Transactions on Cybernetics, 51(6):3115–3128, 2021.
  * [15] Chao Qian.  Distributed pareto optimization for large-scale noisy subset selection.  IEEE Transactions on Evolutionary Computation, 24(4):694–707, Aug 2020.
  * [16] Luis Miguel Antonio and Carlos A. Coello Coello.  Use of cooperative coevolution for solving large scale multiobjective optimization problems.  In 2013 IEEE Congress on Evolutionary Computation, pages 2758–2765, 2013.
  * [17] Cheng He, Lianghao Li, Ye Tian, Xingyi Zhang, Ran Cheng, Yaochu Jin, and Xin Yao.  Accelerating large-scale multiobjective optimization via problem reformulation.  IEEE Transactions on Evolutionary Computation, pages 1–1, 2019\.
  * [18] Hong Qian and Yang Yu.  Solving high-dimensional multi-objective optimization problems with low effective dimensions.  In Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence, AAAI’17, page 875–881. AAAI Press, 2017.
  * [19] Min Jiang, Zhenzhong Wang, Haokai Hong, and Gary G. Yen.  Knee point-based imbalanced transfer learning for dynamic multiobjective optimization.  IEEE Transactions on Evolutionary Computation, 25(1):117–129, 2021\.
  * [20] Ye Tian, Xiutao Zheng, Xingyi Zhang, and Yaochu Jin.  Efficient large-scale multi-objective optimization based on a competitive swarm optimizer.  IEEE Transactions on Cybernetics, pages 1–13, 2019.
  * [21] Min Jiang, Liming Qiu, Zhongqiang Huang, and Gary G. Yen.  Dynamic multi-objective estimation of distribution algorithm based on domain adaptation and nonparametric estimation.  Information Sciences, 435:203 – 223, 2018.
  * [22] M. Jiang, Z. Huang, G. Jiang, M. Shi, and X. Zeng.  Motion generation of multi-legged robot in complex terrains by using estimation of distribution algorithm.  In 2017 IEEE Symposium Series on Computational Intelligence (SSCI), pages 1–6, 2017.
  * [23] X. Ma, F. Liu, Y. Qi, X. Wang, L. Li, L. Jiao, M. Yin, and M. Gong.  A multiobjective evolutionary algorithm based on decision variable analyses for multiobjective optimization problems with large-scale variables.  IEEE Transactions on Evolutionary Computation, 20(2):275–298, 2016\.
  * [24] Wenjing Hong, Ke Tang, Aimin Zhou, Hisao Ishibuchi, and Xin Yao.  A scalable indicator-based evolutionary algorithm for large-scale multiobjective optimization.  IEEE Transactions on Evolutionary Computation, 23(3):525–537, June 2019.
  * [25] Zhenyu Yang, Ke Tang, and Xin Yao.  Large scale evolutionary optimization using cooperative coevolution.  Information ences, 178(15):2985–2999, 2014.
  * [26] H. Zille, H. Ishibuchi, S. Mostaghim, and Y. Nojima.  A framework for large-scale multiobjective optimization based on problem transformation.  IEEE Transactions on Evolutionary Computation, 22(2):260–275, 2018\.
  * [27] Ruochen Liu, Jin Liu, Yifan Li, and Jing Liu.  A random dynamic grouping based weight optimization framework for large-scale multi-objective optimization problems.  Swarm and Evolutionary Computation, 55:100684, 2020.
  * [28] F. van den Bergh and A. P. Engelbrecht.  A cooperative approach to particle swarm optimization.  IEEE Transactions on Evolutionary Computation, 8(3):225–239, 2004\.
  * [29] R. Shang, K. Dai, L. Jiao, and R. Stolkin.  Improved memetic algorithm based on route distance grouping for multiobjective large scale capacitated arc routing problems.  IEEE Transactions on Cybernetics, 46(4):1000–1013, 2016.
  * [30] Minghan Li and Jingxuan Wei.  A cooperative co-evolutionary algorithm for large-scale multi-objective optimization problems.  In Proceedings of the Genetic and Evolutionary Computation Conference Companion, GECCO ’18, page 1716–1721, New York, NY, USA, 2018. Association for Computing Machinery.
  * [31] Y. Yuan, H. Xu, B. Wang, and X. Yao.  A new dominance relation-based evolutionary algorithm for many-objective optimization.  IEEE Transactions on Evolutionary Computation, 20(1):16–37, 2016\.
  * [32] X. Zou, Y. Chen, M. Liu, and L. Kang.  A new evolutionary algorithm for solving many-objective optimization problems.  IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), 38(5):1402–1412, 2008.
  * [33] David Hadka.  Borg: An auto-adaptive many-objective evolutionary computing framework.  Evolutionary Computation, 21(2):231, 2013.
  * [34] F. di Pierro, S. Khu, and D. A. Savic.  An investigation on preference order ranking scheme for multiobjective evolutionary optimization.  IEEE Transactions on Evolutionary Computation, 11(1):17–45, 2007\.
  * [35] Gaoping Wang and Huawei Jiang.  Fuzzy-dominance and its application in evolutionary many objective optimization.  In 2007 International Conference on Computational Intelligence and Security Workshops (CISW 2007), pages 195–198, 2007.
  * [36] Shufen Qin, Chaoli Sun, Yaochu Jin, Ying Tan, and Jonathan Fieldsend.  Large-scale evolutionary multi-objective optimization assisted by directed sampling.  IEEE Transactions on Evolutionary Computation, pages 1–1, 2021\.
  * [37] S. Yang, M. Li, X. Liu, and J. Zheng.  A grid-based evolutionary algorithm for many-objective optimization.  IEEE Transactions on Evolutionary Computation, 17(5):721–736, 2013\.
  * [38] Mario Köppen and Kaori Yoshida.  Substitute distance assignments in nsga-ii for handling many-objective optimization problems.  In Shigeru Obayashi, Kalyanmoy Deb, Carlo Poloni, Tomoyuki Hiroyasu, and Tadahiko Murata, editors, Evolutionary Multi-Criterion Optimization, pages 727–741, Berlin, Heidelberg, 2007. Springer Berlin Heidelberg.
  * [39] J. Cheng, G. G. Yen, and G. Zhang.  A many-objective evolutionary algorithm with enhanced mating and environmental selections.  IEEE Transactions on Evolutionary Computation, 19(4):592–605, 2015\.
  * [40] R. Wang, R. C. Purshouse, and P. J. Fleming.  Preference-inspired coevolutionary algorithms for many-objective optimization.  IEEE Transactions on Evolutionary Computation, 17(4):474–494, 2013\.
  * [41] X. Zhang, Y. Tian, and Y. Jin.  A knee point-driven evolutionary algorithm for many-objective optimization.  IEEE Transactions on Evolutionary Computation, 19(6):761–776, 2015\.
  * [42] C. He, S. Huang, R. Cheng, K. C. Tan, and Y. Jin.  Evolutionary multiobjective optimization driven by generative adversarial networks (gans).  IEEE Transactions on Cybernetics, pages 1–14, 2020.
  * [43] R. Cheng, Y. Jin, M. Olhofer, and B. sendhoff.  Test problems for large-scale multiobjective and many-objective optimization.  IEEE Transactions on Cybernetics, 47(12):4108–4121, 2017.
  * [44] K. Deb and H. Jain.  An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part i: Solving problems with box constraints.  IEEE Transactions on Evolutionary Computation, 18(4):577–601, 2014\.
  * [45] Q. Zhang, A. Zhou, and Y. Jin.  Rm-meda: A regularity model-based multiobjective estimation of distribution algorithm.  IEEE Transactions on Evolutionary Computation, 12(1):41–63, 2008\.
  * [46] Y. Tian, R. Cheng, X. Zhang, and Y. Jin.  Platemo: A matlab platform for evolutionary multi-objective optimization [educational forum].  IEEE Computational Intelligence Magazine, 12(4):73–87, 2017.
  * [47] Goshgar Ismayilov and Haluk Rahmi Topcuoglu.  Dynamic multi-objective workflow scheduling for cloud computing based on evolutionary algorithms.  In 2018 IEEE/ACM International Conference on Utility and Cloud Computing Companion (UCC Companion), pages 103–108, 2018.
  * [48] E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, and V. G. da Fonseca.  Performance assessment of multiobjective optimizers: an analysis and review.  IEEE Transactions on Evolutionary Computation, 7(2):117–132, 2003\.


# The structure of cold neutron star with a quark core within the MIT and NJL
models

T. Yazdizadeh1 111E-mail: tyazdizadeh@yahoo.comr,G. H. Bordbar
2,3,4222Corresponding Author. email<EMAIL_ADDRESS>1Department of
Physics, Payame Noor University (PNU), P.O. Box 19395-3697, Tehran, Iran
2Department of Physics, Shiraz University, Shiraz 71454, Iran333Permanent
address
3Research Institute for Astronomy and Astrophysics of Maragha, PO Box
55134-441, Maragha, Iran
4 Department of Physics and Astronomy, University of Waterloo, 200 University
Avenue West, Waterloo, Ontario N2L3G1, Canada

###### Abstract

Neutron star due to their high interior matter density are expected to be
composed of a quark core, a mixed quark-hadron matter, and a layer of hadronic
matter. Thus, in this paper, we compute the equation of state of these parts
of neutron star to evaluate its structure properties. We use two models for
describing EOS of quark matter, NJL and MIT bag models, and employ three
approaches in this work. A density dependent bag constant satisfy the quark
confinement in the simple MIT bag model. We also study the interaction
behavior of quarks, firstly one gluon exchange within MIT bag model and the
secondly dynamical mass will be held as effective interaction that roles
between particles. Density dependence of quark mass is obtained from NJL self
consistent model. NJL model is a effective manner for justify the chiral
symmetry. Applying the Gibbs conditions the equation of state of the quarks
and hadrons mixed phase is obtained. Since the hadronic matter is under the
influence of strong force of nucleons, we calculate the equation of state of
this phase using a powerful variational many-body technique. Finally, we
calculate the mass and radius of a cold neutron star with a quark core by
numerically solving the TOV equation. To check our used EOS, we compare our
results with the recent observational data. Our results are in a good
agreement with some observed compact objects such as $SAXJ1748.9-2021$,
$4U1608-52$ and $VelaX-1$.

###### pacs:

21.65.-f, 26.60.-c, 64.70.-p

Keywords: neutron star, quark core, quark matter, MIT bag model, NJL model,
gravitational mass, radius.

##  I INTRODUCTION

Neutron stars are placed in category of compact objects, and the interior
matter of them can reach a density much greater than the normal nuclear
saturation density. Therefore, these astrophysical objects are best laboratory
and a unique environment to probe the properties of dense matter. Studying
these stars is one of main problems in physics. In high densities, hadrons
dissolve to quarks, and a phase transition is happen from hadronic matter to
quark matter. Many years ago, the presence of quark matter in neutron stars
has been suggested by Ivanenko1 , Itoh 2  and Collins 3 . There are up, down,
and strange quarks in the quark matter, and this strange matter is a fermi gas
which the other quarks because of their high masses do not appear in this
part. Since all the hadrons do not converge to quarks simultaneously, it is
expected to exist a mixed phase of quarks and hadrons at finite range of
density that the energy is lower than that of quark and hadron matters.

Historically Glendenning was the first who pointed to the neutrality charge of
two phases in mixed phase 4 . In mixed phase, we study the transition from a
hadron phase to quark phase using the Gibbs conditions. Since the existence of
mixed phase of quarks and hadrons affects the properties of neutron star, we
consider the neutron stars to be composed of a quark matter core, a mixed
phase of quarks and hadrons and a layer of hadrons. There is a high
uncertainty in equation of state (EOS) of quark matter. Usually, two more
efficient models are used to study deconfined quark matter, the MIT bag model
5 ; 6  and the Nambu-Jona-Lasinio (NJL) model 7 . The total energy density in
MIT bag model is the sum of the kinetic energy of free quarks and a bag
constant $B$ that is nonperturbative energy shift. We will be search influence
of one-gluon-exchange in MIT bag model. There is asymptotically interaction
among quarks at high densities. This interaction can show up by one gluon
exchange. Therefore, in addition of $B$, we add another term to EOS that is
identified by $\alpha$ the QCD coupling constant 8 ; 9 ; 10 . In MIT bag
model, the mass of quarks is constant, but in the second method, NJL, mass of
quarks depends on density that is considered as the effective interaction of
quarks. Historically first time the NJL model is presented in two papers from
Nambu and Jona-Lasinio in 1961 11 ; 12 . In spite of MIT bag model, the NJL
model does not have confinement, but it justifies chiral symmetry. But at high
density matter like quark matter, although the confinement and chiral symmetry
are in the least importance, both models treat to be similar 13 .

In recent years, we have applied the MIT bag model to investigate the cold and
hot strange quark star. For example we consider MIT bag model with a density
dependent bag constant for a hot strange star in 14 , and found that the mass
and radius of strange star decreases when temperature increases. We have also
found that a higher mass and radius will be obtained for a density dependent
bag constant compare to a fixed bag constant. In 15 ; 16  using MIT bag model
with fixed bag constant, we considered the stability of spin polarized quark
star in a strong magnetic field compared to unpolarized case, and calculated
the structure of this star at zero and finite temperatures. Also We have used
MIT bag model with a density dependent bag constant to calculate the structure
of spin polarized strange star in presence of magnetic field at zero and
finite temperatures 17 ; 18 . We have also applied NJL model to calculate the
equation of state of quark matter 19 .

As we know, a neutron star with a quark core is called hybrid star. In our
previous works, we have obtained the structure of hybrid star by MIT bag model
with fixed and density dependent bag constant model at zero temperature 20
and at finite temperature 21 ; 22 . In those works, we have considered the
simplest version of the MIT bag model. In the present work, we intend to
develop our previous calculations by considering the effects of one gluon
exchange for the quark matter in a neutron star with a quark core. We also use
the NJL method in this work. NJL model is an effective theory and a good
choice for the studying the chiral quark and diquark condensates. Chiral
symmetry and its breakdown in vacuum are the basic property of NJL model. The
outline of our work is as follows: In section II, we calculate the equation of
state of three mentioned phases of the neutron star matter. Then using this
equation of state, we determine the mass and radius of neutron star with a
quark core in section III.

##  II equation of state of a hybrid neutron star

Here, we determine EOS of different part of neutron star: a hadron phase,
quark phase and a mixed phase of quarks and hadrons respectively.

###  II.1 Hadron Phase

We consider the lowest order constrained variational (LOCV) many-body method
for hadron phase 23 ; 24 ; 25 ; 26 . By considering a many body trail wave
function such as $\psi=F\phi$ and some calculations, the cluster expansion of
the energy functional is gotten 31 ,

$\displaystyle E([f])$ $\displaystyle=$
$\displaystyle\frac{1}{A}\frac{<\psi\mid H\mid\ \psi>}{<\psi\mid\psi>}$ (1)
$\displaystyle=$ $\displaystyle E_{1}+E_{2}+E_{3}+\cdots,$

where $F=S\prod_{i<j}f(ij)$ is an A-body correlation operator ($f(ij)$ is two-
body correlation function and $S$ leads to a symmetric product) and $\phi$ is
the slater determinate of $A$ noninteracting nucleons. We consider the first
two terms in above equation, the one-body term
$E_{1}=\sum_{i=1,2}{\frac{3}{5}\frac{k_{i}^{2}}{2m_{i}}\frac{\rho_{i}}{\rho}}$,
and the two-body term $E_{2}=\frac{1}{2A}\sum_{ij}<ij|\nu(12)|ij-ji>$. In
these relations, $\rho_{i}$ is the nucleon density, $\rho=\rho_{p}+\rho_{n}$
is the total density and $k_{i}=(3\pi^{2}\rho_{i})^{1/3}$. Here
$\nu=-\frac{\hbar^{2}}{2m}[f(12),[{\nabla_{12}}^{2},f(12)]]+f(12)V(12)f(12)$
is the nucleonic effective potential ($V(12)$ is the nuclear potential). See
reference 24  for full nuclear matter calculations.

###  II.2 Quark Phase

We employ MIT bag model and NJL model which are two well-known and efficient
models for describing the characteristics of deconfined quark matter.

####  II.2.1 The MIT Bag Model

Quark matter is a fermi gas which composed of deconfined up, down, and strange
quarks. Therefore the total energy is given by

${\cal E}_{tot}={\cal E}_{u}+{\cal E}_{d}+{\cal E}_{s}+B.$ (2)

In equation (2) ${\cal E}_{i}$ is

$\displaystyle{\cal E}_{i}$ $\displaystyle=$
$\displaystyle\frac{3m_{i}^{4}}{8\pi^{2}}\left[x_{i}(2x_{i}^{2}+1)(\sqrt{1+x_{i}^{2}})-\sinh^{-1}x_{i}\right]$
(3)
$\displaystyle-\alpha_{c}\frac{m_{i}^{4}}{\pi^{3}}\left[x_{i}^{4}-\frac{3}{2}[x_{i}(\sqrt{1+x_{i}^{2}})-\sinh^{-1}x_{i}]^{2}\right],$

where

$x_{i}=\frac{k_{F}^{(i)}}{m_{i}}.$ (4)

Here, $k_{F}^{(i)}=(\rho_{i}\pi^{2})^{1/3}$ , $m_{i}$ and $\rho_{i}$ are the
mass and baryon density of quark $i$, respectively, and $\alpha_{c}$ is the
QCD coupling constant. We consider three values for $\alpha$: $\alpha=0$ (MIT
bag model without interaction), $\alpha=0.16$ and $\alpha=0.5$. Although these
amount of $\alpha$ are small and perturbative, they show an appropriate range
of quark interaction and are in the selection rang of Farhi and Jaffe work f .
In equation (2), $B$ is a density dependent bag constant which satisfies the
quark confinement in MIT bag model. We consider a Gaussian form
$B(\rho)=B_{\infty}+(B_{0}-B_{\infty})\exp\left[{-\beta(\frac{\rho}{\rho_{0}})^{2}}\right]$,
where $B_{0}=B(\rho=0)=400MeV/fm^{3}$, $\beta$ is a numerical parameter equal
to $\rho_{0}=0.17fm^{-3}$ and $B_{\infty}$ is a free parameter which is
determined by using the experimental data reported in the CERN SPS. 20 ; 34 .
Now, using the energy density from Eq. (2), the EOS of quark matter in the MIT
bag model is obtained,

$P(\rho)=\rho\frac{\partial{\cal E}}{\partial\rho}-{\cal E}.$ (5)

####  II.2.2 The NJL Model

In NJL model, The dynamical mass is held the effective interaction between
particles. In this method, we adopt a lagrangian similar to that given in
reference 19 , as follow,

$\displaystyle{\cal L}$ $\displaystyle=$
$\displaystyle\bar{q}\left(i\gamma^{\mu}\partial_{\mu}-\hat{m}_{0}\right)q+G\sum_{k=0}^{8}\left[(\bar{q}\lambda_{k}q)^{2}+(\bar{q}i\gamma_{5}\lambda_{k}q)^{2}\right]$
(6)
$\displaystyle-K\left[det_{f}(\bar{q}(1+\gamma_{5})q)+det_{f}(\bar{q}(1-\gamma_{5})q)\right],$

where $q$ is the field of quarks with three flavors and three colors,
$\hat{m}_{0}=diag(m_{0}^{u},m_{0}^{d},m_{0}^{s})$ in flavor space, and
$\lambda_{k}$ $(0\leq k\leq 8)$ are the flavor matrices. Restoring chiral
symmetry, breaking is indicated with a ultra-violet cut-off. We employ
parameters of reference 32 ; 33  as follow $\Lambda=602.3MeV$. ;
$G\Lambda^{2}=1.835\Lambda$ and $K\Lambda^{2}=12.36$. $G$ and $K$ are coupling
strength.

The dynamical mass is calculated by

$m_{i}=m_{0}^{i}-4G<\bar{q}_{i}q_{i}>+2K<\bar{q}_{j}q_{j}><\bar{q}_{k}q_{k}>,$
(7)
$<\bar{q}_{i}q_{i}>=-\frac{3}{\pi^{2}}\int_{p_{Fi}}^{\Lambda}p^{2}dp\frac{m_{i}^{2}}{\sqrt{m_{i}^{2}+p^{2}}},$
(8)

$p_{Fi}=(\pi^{2}\rho_{i})^{1/3}.$ (9)

At this stage we will determine EOS of quark matter in NJL model,

$P=\sum_{i=u,d,s}n_{i}\sqrt{p_{Fi}^{2}+m_{i}^{2}}-{\cal E},$ (10)

where

${\cal
E}=\sum_{i=u,d,s}\frac{3}{\pi^{2}}\int_{0}^{p_{Fi}}p^{2}dp\sqrt{p^{2}+m_{i}^{2}}-(B-B_{0}).$
(11)

In equation (11), $B$ is the bag pressure 13  which is consequence of
interaction,

$\displaystyle B$ $\displaystyle=$
$\displaystyle\sum_{i=u,d,s}\left[\frac{3}{\pi^{2}}\int_{0}^{\Lambda}p^{2}dp\left(\sqrt{p^{2}+m_{i}^{2}}-\sqrt{p^{2}+{m_{0}^{i}}^{2}}\right)\right.$
(12)
$\displaystyle\left.-2G<\bar{q}_{i}q_{i}>^{2}\right]+4K<\bar{u}u><\bar{d}d><\bar{s}s>.$

###  II.3 Mixed phase

There is a mixed phase of quarks and hadrons within a finite rang of density.
In this phase, we apply the Gibss condition. According to this equilibrium
condition, the pressures and chemical potentials of both quark and hadron
phases are equal 4 ,

$\mu^{Q}_{n}=\mu^{H}_{n},$ (13) $\mu^{Q}_{p}=\mu^{H}_{p},$ (14)

where $\mu^{H}_{n}$ ($\mu^{H}_{p}$) and $\mu^{Q}_{n}$ ($\mu^{Q}_{p}$) are the
neutrons (protons) chemical potential in nucleonic and quark part in the mixed
phase, respectively. Using above equations, we can determine the charge
density of quarks and hadrons. $\chi$, which is the volume fraction occupied
by quarks is determined by considering the global charge neutrality. Then
baryon and total energy densities in mixed phase could be calculated.

$\chi(\frac{2}{3}\rho_{u}-\frac{1}{3}\rho_{d}-\frac{1}{3}\rho_{s})+(1-\chi)\rho_{p}-\rho_{e}=0,$
(15) $\rho_{B}=\chi\rho_{Q}+(1-\chi)\rho_{H},$ (16) ${\cal E}_{MP}=\chi{\cal
E}_{QP}+(1-\chi){\cal E}_{HP}.$ (17)

See reference 41  for detail of calculations regarding the EOS of mixed phase.

At this stage we can determine EOS of neutron star with a quark core using the
results of proceeding sections. We consider three approaches for the quark
matter (in the quark phase and mixed phase); simple MIT bag model (model 1),
MIT bag model by considering the quark interaction (model 2) and NJL model
(model 3). EOS results for neutron star (corresponding to different models for
the quark matter calculations) are given in Fig. 1.

Figure 1: EOS for the hybrid neutron star. The results of our calculations in
different models of quark matter have been plotted.

This figure shows that the difference between the equation of state in model 1
and model 3 is substantially lower than those of other mentioned models.
However, we see that the equation of state of neutron star becomes stiffer
when coupling constant is higher, particulary at high densities. This shows
the importance of interaction at higher densities.

##  III Structure properties of the hybrid Neutron Star

Using the obtained EOS in previous section, the structure of neutron star with
a quark core can be calculated. Before this work, we investigate the energy
and stability conditions for our results. For this purpose at first, we fit a
polynomial function for the equations of state in Fig. 1 as ,

$P=\sum_{i=1}^{7}a_{i}{\cal E}^{7-i}.$ (18)

The coefficients $a_{i}$ have been given in Table 1.

Table 1: Different coefficients presented in Eq. (18) for our applied models. Model | $a_{1}\ (\times 10^{-57})$ | $a_{2}\ (\times 10^{-40})$ | $a_{3}\ (\times 10^{-25})$ | $a_{4}\ (\times 10^{-10})$ | $a_{5}\ (\times 10^{6})$ | $a_{6}\ (\times 10^{21})$ | $a_{7}\ (\times 10^{35})$
---|---|---|---|---|---|---|---
model 1 | $1.194$ | $-0.2467$ | $2.011$ | $-8.123$ | $1.656$ | $-1.201$ | $2.915$
model 2; ${\alpha=0.16}$ | $6.487$ | $-1.279$ | $9.852$ | $-37.23$ | $7.090$ | $-4.962$ | $10.80$
model 2; ${\alpha=0.5}$ | $9.400$ | $-1.912$ | $15.11$ | $-58.19$ | $11.20$ | $-6.819$ | $11.27$
model 3 | $1.285$ | $-0.2407$ | $1.793$ | $-6.630$ | $1.209$ | $-0.7027$ | $1.321$

We use this relation to justify the energy and stability conditions as
follows.

###  III.1 Energy conditions

Energy conditions in the center of neutron star include the null energy
condition (NEC), weak energy condition (WEC), strong energy condition (SEC)
and dominant energy condition (DEC). These condition are expressed as follow,

$NEC\rightarrow P_{c}+{\cal E}_{c}\geq 0,$ (19)

$WEC\rightarrow P_{c}+{\cal E}_{c}\geq 0\quad\&\quad{\cal E}_{c}\geq 0,$ (20)

$SEC\rightarrow P_{c}+{\cal E}_{c}\geq 0\quad\&\quad 3P_{c}+{\cal E}_{c}\geq
0,$ (21)

$DEC\rightarrow{\cal E}_{c}>|P_{c}|,$ (22)

where ${\cal E}_{c}$ is the energy density and $P_{c}$ is the pressure at the
center of star. Results of the above conditions for our equations of state are
given in Table 2.

Table 2: Energy conditions for hybrid neutron star for applied models. | NS+Quark Core | ${\cal E}_{c}\left({10^{14}gr/cm^{3}}\right)$ | $P_{c}\left({10^{14}gr/cm^{3}}\right)$ | $NEC$ | $WEC$ | $SEC$ | $DEC$
---|---|---|---|---|---|---|---
| model 1 | 25.8 | 7.79 | $\surd$ | $\surd$ | $\surd$ | $\surd$
| model 2 ; ${\alpha=0.16}$ | 18.25 | 17.8 | $\surd$ | $\surd$ | $\surd$ | $\surd$
| model 2 ; ${\alpha=0.5}$ | 14.05 | 28.9 | $\surd$ | $\surd$ | $\surd$ | $\times$
| model 3 | 22.9 | 5.56 | $\surd$ | $\surd$ | $\surd$ | $\surd$

We have found that our equations of state satisfy mentioned energy conditions,
except the dominant energy condition for model 2 with $\alpha=0.5$ in which
the central pressure is very high with respect to the other models.

###  III.2 Stability

According to the stability condition, an equation of state is physically
acceptable when the corresponding obtained velocity of sound ($v$) be less
than the light’s velocity (c) 35 ; 36 . Thus the stability condition is
$(0\leq v^{2}=(\frac{dP}{d\cal E})\leq c^{2})$. By Using Eq. 18, we have
computed $\frac{v^{2}}{c^{2}}$ versus density which has been given in Fig. 2.
It is evident that the stability condition is satisfied by the our calculated
EOS of neutron star with quark core for models 1 and 3.

Figure 2: Sound speed versus density for different models.

So these two models are suitable for determining the structure of a hybrid
neutron star. The other models do not obey the stability condition, and we
can’t use their equations of state in the structure calculations.

###  III.3 Properties of the neutron star with a quark core

We use TOV equation to calculate the structure of star. This equation is
determined by Tolman-Oppenheimer-Volkoff (TOV) 37 ; 38 ; 39 .

$\frac{dP}{dr}=-\frac{G\left[{\frac{P}{c^{2}}+\cal E}\right]\left[m+\frac{4\pi
r^{3}P}{c^{2}}\right]}{r^{2}\left[1-\frac{2Gm}{rc^{2}}\right]},$ (23)

$\frac{dm}{dr}=4\pi r^{2}{\cal E}$ (24)

where $P$ is the pressure and $\cal E$ is the energy density. In our
calculations for the neutron star, we consider the following equation of
state: up to the density $0.05fm^{-3}$, we use the data of Baym calculations
40 . For the hadron, quark, and mixed phases, we consider our equations of
state obtained in the previous sections. Numerically integrating the TOV
equation for a given equation of state, the mass and radius of the neutron
star with a quark core is determined. Our results are as follows.

The gravitational mass of hybrid neutron star (a neutron star with a quark
core) have been presented in Fig. 3 versus the central mass density. The mass-
radius relation has been also shown in Fig. 4 for this star for different
models. For comparison, we have also brought the results for a neutron star
without a quark matter 41  in Figs. 3 and 4.

Figure 3: Gravitational mass versus the central mass density for the neutron
star with a quark core for different models. The results for a neutron star
without quark core (NS) have been also given for comparison. Figure 4: Mass-
radius relation for the neutron star with a quark core for different models.
The results of neutron star without quark core (NS) have been also given for
comparison.

Results show that in the NJL model with a density dependent mass of quarks,
the maximum mass is lower than that of other models. The structure properties
of neutron star in the cases without and with a quark and have been given in
Table 3 for different models.

Table 3: Structure properties of neutron star without (NS) and with (NS+Q) a quark core for different models. The observational data have been also presented for comparison. | Compact object | $M_{max}\left(M_{\,\odot}\right)$ | $R\left(km\right)$
---|---|---|---
| NS | 1.98 | 9.8
| NS+Q: model 1 | 1.8 | 10
| NS+Q: model 3 | 1.75 | 10.4
|  $4U1820-30$ 42  | $1.58\pm 0.6$ | $9.1\pm 0.4$
|  $PSRJ1903+0327$ 43  | $1.667\pm 0.021$ | $9.438km$
|  $PSRJ1614-2230$ 44  | $1.97\pm 0.04$ | $13\pm 2km$
|  $SAXJ1748.9-2021$45  | $1.78\pm 0.3$ | $8.18\pm 1.62km$
|  $4U1608-52$ 46  | $1.74\pm 0.14$ | $9.3\pm 1km$
|  $VelaX-1$ 46  | $1.77\pm 0.08$ | $9.56km$
Table 4: Structure properties of neutron star with a quark core is calculated by other authors. Our results have been also given for comparison. | $Reference$ | $Quark\,Phase$ | $Hadron\,Phase$ | $Mixed\,Phase$ | $used$ | $used$ | $M_{max}\left(M_{\,\odot}\right)$ | $R\left(km\right)$
---|---|---|---|---|---|---|---|---
|  | $model$ | $model$ | $condition$ | $parameter$ | $parameter$ |  |
| present wok | MIT Bag Model | LOCV | Gibbs |  |  | 1.8 | 10
| present work | NJL Model | LOCV | Gibbs |  |  | 1.75 | 10.4
|  |  |  |  | $x_{\sigma}=0$ | $G_{2}=0.006$ | 1.44 | 9.54
|  |  |  |  | $x_{\sigma}=0$ | $G_{2}=0.012$ | 1.89 | 12.55
|  |  |  |  | $x_{\sigma}=0$ | $G_{2}=0.016$ | 2.05 | 12.66
| 46  | Field correlator | Relativistic |  | $x_{\sigma}=0.8$ | $G_{2}=0.006$ | 1.44 | 9.52
|  | Method | Mean Field | Gibbs | $x_{\sigma}=0.8$ | $G_{2}=0.012$ | 1.89 | 12.53
|  |  | Model |  | $x_{\sigma}=0.8$ | $G_{2}=0.016$ | 2.04 | 12.40
|  |  |  |  | $x_{\sigma}=0.6$ | $G_{2}=0.006$ | 1.43 | 9.51
|  |  |  |  | $x_{\sigma}=0.6$ | $G_{2}=0.012$ | 1.73 | 12.00
|  |  |  |  | $x_{\sigma}=0.6$ | $G_{2}=0.016$ | 1.8 | 11.49
|  |  |  |  | $\sigma=0$ | $G_{v}=0$ | 1.91 | 13.09
|  |  |  |  | $\sigma=0$ | $G_{v}=0.2G_{s}$ | 2.05 | 13.00
|  |  |  |  | $\sigma=0$ | $G_{v}=0.4G_{s}$ | 2.13 | 12.77
| 47  | NJL Model | Relativistic |  | $\sigma=10$ | $G_{v}=0$ | 1.94 | 13.3
|  |  | Mean Field | Gibbs | $\sigma=10$ | $G_{v}=0.2G_{s}$ | 2.08 | 13.01
|  |  | Model |  | $\sigma=10$ | $G_{v}=0.4G_{s}$ | 2.15 | 12.77
|  |  |  |  | $\sigma=40$ | $G_{v}=0$ | 2.00 | 13.37
|  |  |  |  | $\sigma=40$ | $G_{v}=0.2G_{s}$ | 2.11 | 13.03
|  |  |  |  | $\sigma=40$ | $G_{v}=0.4G_{s}$ | 2.17 | 12.67
|  |  |  |  |  | $\beta=0.000$ | 2.029 | 11.13
|  |  |  |  |  | $\beta=0.025$ | 2.003 | 11.56
| 48  | MIT | Chiral Quark- | Gibbs | ESO8 | $\beta=0.050$ | 1.958 | 11.71
|  | Bag Model | Meson Coupling |  |  | $\beta=0.100$ | 1.896 | 11.63
|  |  | Model |  |  | $\beta=0.150$ | 1.866 | 11.5
|  |  |  |  |  | $\beta=0.200$ | 1.853 | 11.41

In table 3, it is seen that when we consider a quark core for the neutron star
in both considered models, the maximum mass decreases and radius increases.
This is because the equation of state becomes softer when we consider a quark
matter in the core of neutron star. We have obtained the lowest amount of mass
in NJL model $(mass=1.75M_{\,\odot})$, but it is nearly close to that of MIT
bag model with zero coupling constant $(mass=1.8M_{\,\odot})$. The calculated
radius in NJL model is bigger than that of MIT bag model. We have also brought
some observational data of neutron star candidates in table 3 to compare our
results with these data. The mass of a neutron star with a quark core which we
have been determined in this paper is closed to the mass of $SAXJ1748.9-2021$
46 , $4U1608-52$ 47 , and $VelaX-1$ 47 . Also the calculated radius of neutron
star with a quark core in two models are in a good agreement with the radius
of the mentioned observed compact objects. However, the mass and radius do not
agree with the recent observational data of for pulsar $PSRJ1614-2230$ 45
with $M=1.97\pm 0.04M_{\,\odot}$, while our result of neutron star without
quark core 41  has a good agreement with the mass of this pulsar. But the
calculated radius is smaller than the radius of pulsar $PSRJ1614-2230$.

Also we review the works of several authors who have researched on the
properties of hybrid neutron star. The results of their models are presented
briefly in table 4. In this table, $x_{\sigma}$ is the hyperon coupling
constant and $G_{2}$ is one of the EOS parameters within the field correlated
method 47 . We can see that the results of reference 47  with
$x_{\sigma}=0.6$, and $G_{2}=0.012$ and $0.016$ ($M=1.73M_{\,\odot}$ and
$1.8M_{\,\odot}$) are in a good agreement with our calculated mass. In table
4, $\sigma$ and $G_{v}$ are the surface tension and the vector coupling
constant, respectively 48 . Mass values in that paper are about
$2M_{\,\odot}$, therefore they have a good agreement with $PSRJ1614-2230$,
while they do not have agreement with our work. In table 4, ESCO8 is a model
for determining hyperon coupling constant and $\beta$ is a parameter which is
the same with the one in bag constant formula, $B$, in our work 49 . The
obtained masses in this reference is higher than our calculated result. The
authors in all these papers have parameterized the EOS in order to get $M\cong
2.0M_{\,\odot}$ which has been recently observed for $PSRJ1614-2230$. Here, we
can conclude that if we want to get the mass of a hybrid star to be about the
mass of pulsar $PSRJ1614-2230$ with the mass $M=1.97\pm 0.04M_{\,\odot}$, we
should modify the equation of state that used in this work, or we should use
another computational method to calculate the structure of a hybrid star.

##  IV Summary and Conclusion

Since the neutron stars are one of the compact objects with high density, this
idea raises that there is a deconfined quark matter in these stars. Thus, here
we considered a crust of hadronic matter, a mixed phase of quark and hadronic
matters and a quark matter in the core of neutron star. In this work, we
calculated the EOS for quark matter phase of neutron star in three models,
simple MIT bag model, MIT bag model including one gluon exchange correction
with two different coupling constants and NJL model with a density dependent
mass for quarks. For hadronic matter phase, we chose a variational method
(LOCV). For the mixed phase, the hadron-quark phase transition was modeled by
the Gibbs constructions. After calculation of EOS, we studied the energy and
stability conditions. We found that when the interaction of quarks is given by
the one gluon exchange, the equation of state doesn’t satisfy these
conditions. Therefor in this case EOS isn’t suitable for calculation of
structure of this star. Using determined equation of state and solving Tolman-
Oppenheime-Volkof (TOV) equations, we computed the structure of a neutron star
with a quark core (hybrid star) in two models 1 and 3. We saw that our result
for the maximum mass of neutron star with a quark core agrees with the
observed mass for $SAXJ1748.9-2021$, $4U1608-52$ and $VelaX-1$. However, that
is not in a good agreement with the recent observational data for
$PSRJ1614-2230$. While, we found that our result for the maximum mass of
neutron star without quark matter is in agreement to the mass of this object.
Here, it can be concluded that our equation of state should be modified in
order to get some good agreements with the new observational data.

## Acknowledgements

We wish to thank the Research Council of Islamic Azad University, Bafgh
Branch. We also wish to thank Shiraz University Research Council. This work
has been supported financially by the Research Institute for Astronomy and
Astrophysics of Maragha, Iran. G. H. Bordbar wishes to thank N. Afshordi
(University of Waterloo) for his useful comments and discussions during this
work. G. H. Bordbar also wishes to thank Physics Department of University of
Waterloo for the great hospitality during his sabbatical.

## References

  * (1) D. D. Ivanenko and D. F. Kurdgelaidze, _Astrophys._ 1 , 251 (1965).
  * (2) N. Itoh, _Progr. Theor. Phys._ 44, 291 (1970).
  * (3) J. C. Collins and M. J. Perry, _Phys. Rev. Lett._ 34, 1353 (1975).
  * (4) N. K. Glendenning, _Phys. Rev._ D46, 1274 (1992).
  * (5) A. Chodos, R. L. Jaffe, K. Johnson, C. B. Thorn and V. F. Weisskopf, _Phys. Rev._ D9, 3471 (1974).
  * (6) G. Baym, E. W. Kolb, L. McLerran, T. P. Walker and R. L. Jaffe, _Phys. Lett._ B160, 181 (1985).
  * (7) J. D. Carroll, D. B. Leinweber, A. W. Thomas and A. G. Williams, _Phys. Rev._ C79, 045810 (2009).
  * (8) M. S. Berger and R. L. Jaffe, _Phys. Rev._ C35, 213 (1987).
  * (9) E. Farhi and R. L. Jaffe, _Phys. Rev._ D30, 2379 (1984).
  * (10) T. Maruyama, S. Chiba, H-J. Schulze and T. Tatsumi, _Phys. Rev._ D76, 123015 (2007).
  * (11) Y. Nambu and G. Jona-Lasinio, _Phys. Rev._ 122, 345 (1961).
  * (12) Y. Nambu and G. Jona-Lasinio, _Phys. Rev._ 124, 246 (1961).
  * (13) M. Buballa, _Phys. Rep._ 407, 205 (2005).
  * (14) G. H. Bordbar, A. Poostforush and A. Zamani, _Astrophys._ 54, 277 (2011).
  * (15) G. H. Bordbar and A. Peyvand, _Res. Astron. Astrophys._ 11, 851 (2011).
  * (16) G. H. Bordbar, F. Kayanikhoo and H. Bahri, _Iranian J. Sci. Tech._ A37, 165 (2013).
  * (17) G. H. Bordbar, H. Bahri and F. Kayanikhoo, _Res. Astron. Astrophys._ 12, 1280 (2012).
  * (18) G. H. Bordbar and Z. Alizade, _Astrophys._ 57, 130 (2014).
  * (19) G. H. Bordbar and B. Ziaei, _Res. Astron. Astrophys._ 12, 540 (2012).
  * (20) G. H. Bordbar, M. Bigdeli and T. Yazdizadeh, _Int. J. Mod. Phys._ A21, 5991 (2006).
  * (21) T. Yazdizadeh and G. H. Bordbar, _Res. Astron. Astrophys._ 11, 471 (2011).
  * (22) T. Yazdizadeh and G. H. Bordbar, _Astrophys._ 56, 121 (2013).
  * (23) G. H. Bordbar and M. Modarres, _J. Phys. G: Nucl. Phys._ 23, 1631 (1997).
  * (24) G. H. Bordbar and M. Modarres, _Phys. Rev._ C57, 714 (1998).
  * (25) M. Modarres and G. H. Bordbar, _Phys. Rev._ C58, 2781 (1998).
  * (26) G. H. Bordbar and M. Bigdeli, _Phys. Rev._ C75, 045804 (2007).
  * (27) Fahri, E., and Jaffe, R. L., Phys. Rev., D30, 1601 (1984).
  * (28) J. W. Clark, _Prog. Part. Nucl. Phys._ 2, 89 (1979).
  * (29) P. Rehberg, S. P. Klevansky and J. Haufner, _Phys. Rev._ C53, 410 (1996).
  * (30) M. Buballa and M. oertel, _Phys. Lett._ B457, 261 (1999)
  * (31) G. F. Burgio, M. Baldo, P. K. Sahu and H. J. Schulze, _Phys. Rev._ C66, 025802 (2002).
  * (32) L. Herrera, _Phys. Lett._ A165, 206 (1992).
  * (33) H. Abreu, H. Hernandez and L. A. Nunez, _Class. Quant. Grav._ 24 4631 (2007).
  * (34) S. L. Shapiro and S. A. Teukolski, _Black Holes, White Dwarfs and Neutron Stars_ (New York, 1983).
  * (35) N. K. Glendenning, _Compact Star, Nuclear Physics, Particle Physics, and General Relativity_ (Springer, New York, 2000).
  * (36) F. Weber, _Pulsars as Astrophysical Laboratories for Nuclear and Particle Physics_ (Institute of Physics, Bristol, 1999).
  * (37) G. Baym, C. Pethick and P. Sutherland, _Astrophys. J._ 170, 299 (1971).
  * (38) G. H. Bordbar and M. Hayati, _Int. J. Mod. Phys._ A21, 1555 (2006).
  * (39) T. Guver, P. Wroblewski, L. Camarota, and F. Ozel, _Astrophys. J._ 719,1807 (2010).
  * (40) P. C. C. Freire et al., _Mont. Not. R. Astron. Soc._ 412, 2763 (2011).
  * (41) P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J.W. T. Hessels, _Nature (London)_ 467, 1081 (2010).
  * (42) T. Guver and F. Ozel, _Astrophys. J._ 765, L1 (2013)
  * (43) M. L. Rawls et al, _Astrophys. J._ 730, 25 (2011).
  * (44) D. logoteta and I. Bambaci, _Phys. Rev._ D88, 063001 (2013).
  * (45) X. H. Wu and H. Shen, arxiv:1708.01878
  * (46) T. Miyatsu, M. Cheoun and K. Saito, _Astrophys. J._ 813, 135 (2015).