# A Multicriteria Macroeconomic Model with
Intertemporal Equity and Spatial Spillovers

Herb Kunze111Department of Mathematics and Statistics, University of Guelph,
Guelph, Canada. Email<EMAIL_ADDRESS>Davide La Torre222SKEMA Business
School - Universit$\acute{e}$ C$\hat{o}$te d’Azur, Sophia Antipolis Campus,
France. Email<EMAIL_ADDRESS>Simone Marsiglio333Department of
Economics and Management, University of Pisa, Pisa, Italy. Email:
<EMAIL_ADDRESS>
###### Abstract

We analyze a macroeconomic model with intergenerational equity considerations
and spatial spillovers, which gives rise to a multicriteria optimization
problem. Intergenerational equity requires to add in the definition of social
welfare a long run sustainability criterion to the traditional discounted
utilitarian criterion. The spatial structure allows for the possibility of
heterogeneiity and spatial diffusion implies that all locations within the
spatial domain are interconnected via spatial spillovers. We rely on different
techniques (scalarization, $\epsilon$-constraint method and goal programming)
to analyze such a spatial multicriteria problem, relying on numerical
approaches to illustrate the nature of the trade-off between the discounted
utilitarian and the sustainability criteria.

Keywords: Multicriteria Optimization; Intergenerational Equity; Spatial
Diffusion

##  1 Introduction

After decades of debates a wide consensus on the effects of anthropogenetic
activities on the environment has finally emerged, and even policymakers seem
finally convinced that it is now time to act in order to ensure the long term
sustainability of economic activities. Sustainability is a complicated notion
to define but in its most widely accepted terms it requires some respect of
natural resources and some efforts to ensure intergenerational equity ([30]).
While the former aspect can be easily accounted for in macroeconomic analysis
by adding an additional constraint to the standard optimization problem
representing the society’s planning device, the latter is more problematic
since the criterion generally used in the definition of such a problem, the
discounted utilitarianism, attaches less and less weight to future generations
and so it cannot accommodate for intergenerational issues ([13]). Different
approaches to overcome such an issue have been proposed, but probably the most
commonly used consists of extending the objective function to include an
additional term representing somehow long term sustainablity considerations
([12]). The introduction of such an additional term in the optimization
problem makes it a multiple objective problem in which the society needs to
balance two conflicting goals, represented by short term and long term
objectives, respectively; such an interpretation of a macroeconomic problem as
a multicriteria problem allows to rely upon the operations research methods to
analyze macroeconomic issues, bridging somehow these two important but
distinct fields of research ([21]; [22]).

The goal of this paper is therefore to analyze from an operations research
perspective a traditional macroeconomic model extended along two different
directions. The first consists of allowing for long term sustainability
considerations, in the form of intertemporal equity issues, which introduces a
second criterion in the otherwise standard unicriterion macroeconomic
optimization problem; several studies have analyzed in different contexts
problems of this kind by focusing on pollution control and debt reduction
settings ([17]; [19]). The second consists of allowing for a spatial dimension
in which different locations within the entire economy interact one another
through the trade channel; several studies analyze problems of this kind,
which generally are all unicriterion, by focusing on capital accumulation and
environmental problems ([16]; [18]). Specifically, we consider a simple
setting with intergenerational equity (as in [13]) in which capital evolves
both over time and across space (as in [3]), and to the best of our knowledge
our paper represents the first attempt to analyze a spatial macroeconomic
problem from an operational research point of view.

Our paper relates thus to two very distinct literatures addressing the
intertemporal equity issues associated with sustainabiility and the existence
of spatial spillovers in macroeconomic geography contexts, respectively.
Several papers discuss the problems embedded in the use of discounted
utilitarianism as a social welfare criterion by proposing alternative criteria
and discussing their limits in terms of applicability and in terms of
existence of optimal solution paths ([26]; [11]; [23]; [20]; [2]; [14]).
Several other papers instead more recently discuss how the presence of a
spatial dimension allows to characterize the possibility of spatial
heterogeneity and spatial spillovers, along with their implications for
macroeconomic outcomes ([6]; [3]; [7]; [8]; [16]). We bridge these two
branches of the economics literature by developing a spatial version of the
simplest intertemporal equity problem to show how an operations research
approach can be used to inform our analysis.

The paper proceeds as follows. Section 2 recalls the basic definitions and
properties of optimization theory with multiple objectives, discussing the
main methods that we shall adopt in our analysis, that are scalarization,
$\epsilon$-constraint method and goal programming. Section 3 presents our
multicriteria problem showing how it can be reformulated from the points of
view of scalarization, $\epsilon$-constraint method and goal programming.
Section 4 determines the solution of the problem relying on scalarization
illustrating the nature of the trade-off between the discounted utilitarian
and the sustainability criteria and deriving the Pareto-frontier in a specific
model’s parametrization. Section 5 presents some further numerical experiments
by analyzing the problem via the $\epsilon$-constraint method and goal
programming, under different model’s parametrizations. Section 6 as usual
concludes.

##  2 Preliminaries on Multiple Objective Optimization

This section recalls some basic facts in Multiple Objective Optimization
(MOP). In a very abstract setting, a finite-dimensional MOP problem (see [27])
takes the form

$\displaystyle\max_{x\in X}\ J(x)$ (1)

where $(X,\|\cdot\|)$ is a Banach space and $J:X\to{\mathbb{R}}^{p}$ is a
vector-valued functional. As usual we suppose that ${\mathbb{R}}^{p}$ is
ordered by the classical Pareto cone ${\mathbb{R}}^{p}_{+}$. A point $x\in X$
is said to be Pareto optimal or efficient if $J(x)$ is one of the maximal
elements of the set of achievable values $J(X)$. Thus a point $x$ is Pareto
optimal if it is feasible and, for any possible $x^{\prime}\in X$,
$J(x)\leq_{{\mathbb{R}}^{p}_{+}}J(x^{\prime})$ implies $x=x^{\prime}$. In a
more synthetic way, a point $x\in X$ is said to be Pareto optimal if
$(J(x)+{\mathbb{R}}^{p}_{+})\cap J(X)=\\{J(x)\\}$. Because of its
dimensionality and the existing of conflicting criteria, a MOP model is
usually difficult to be solved and the determination of the entire or part of
the Pareto frontier can be very complicated and computationally intensive. In
particular this applies when the number of objectives is larger than two,
leading to a higher-dimensional Pareto surface. To overcome this difficulties
and reduce the model complexity, several techniques have been proposed in
literature. The generation of the Pareto frontier can be accomplished through
one of two predominant techniques: scalarization and vectorization methods.
Among the scalarization techniques, the most frequently applied are linear
scalarization, the $\epsilon$-constraint method, and goal programming. These
techniques will be used in the sequel of this paper. The vectorization
algorithms, instead, tackle the MOP model directly without transforming it
into some equivalent single criterion model.

###  2.1 Linear Scalarization

The linear scalarization technique (or weighted sum) is probably the simplest
and the most widely used technique to solve MOP problems and it converts the
MOP model into a family of parametric single criterion optimization models. By
using this approach, a multiple objective model can be reduced to a single
criterion problem by summing up all criteria with different weights. More
precisely, by linear scalarization a MOP model boils down to:

$\displaystyle\max_{x\in X}\sum_{i=1}^{p}\beta_{i}J_{i}(x)$ (2)

where $\beta$ is a vector taking values in the interior of
${\mathbb{R}}^{p}_{+}$, namely $\beta\in{\rm int}({\mathbb{R}}^{p}_{+})$.
Since the Pareto optimal solution depends on $\beta$, by modifying the weights
$\beta$ different points on the Pareto optimal set can be found. Linear
scalarization can also be applied to problems in which the ordering cone is
different from the Pareto one. In this case, we have to rely on the elements
of the dual cone to scalarize the problem.

If $J$ is a vector-valued concave functional, namely each component $J_{i}$ is
concave, then the linear scalarized problem (2) is also concave. This means
that one can find Pareto optimal points of a concave problem by solving a
concave scalar optimization problem, and for each $\beta\in{\rm
int}({\mathbb{R}}^{p}_{+})$, different Pareto optimal points can be obtained.
For concave problems the converse of this result is only partially true, since
for each Pareto optimal point $\bar{x}$, there is a nonzero
$\bar{\beta}\in{\mathbb{R}}^{p}_{+}$ such that $\bar{x}$ is a solution of the
scalarized problem (2) with $\beta=\bar{\beta}$. This is stated in the
following theorem.

######  Proposition 1.

([27]) Suppose that $D$ is convex and $J_{i}$ are concave for all $i=1...p$.
Then for all Pareto optimal solutions $\bar{x}$ there exists
$\beta\in{\mathbb{R}}^{p}_{+}$ such that

$\bar{x}\in argmax_{x\in D}\left\\{\sum_{i=1}^{p}\beta_{i}J_{i}(x)\right\\}.$
(3)

If the linear scalarization method is used for non-concave problems, the
Pareto frontier generated will be incomplete and the Decision Maker (DM) will
have a non-complete set of possible solutions. In this case, other
scalarization methods can be found in literature and one which is worth to be
mentioned is the Chebyshev scalarization model.

####  2.1.1 $\epsilon$-constraint method

The second model that is proposed to solve the vector-valued problem is the
$\epsilon$-constraint method. In this method, one of the objective functions
is selected to be optimized while the others are converted into additional
constraints. The method is an hybrid methodology, in fact for the
$\\{J_{i}\\}_{i\neq k}$, least acceptable levels, $\epsilon_{i}$ have to be
set while the remaining objective function $J_{k}$ is optimized. Then the
decision maker plays a crucial role in this setting, by choosing which
objective function to optimize and the least acceptable levels for the
objective functions to add as constraints. Under this method, the original
vector-valued problem can be now written as:

$\max J_{k}(x)$ (4)

subject to:

$\left\\{\begin{array}[]{ll}J_{i}(x)\geq\epsilon_{i}\ \ \ \ i\neq k\\\ x\in
X\end{array}\right.$ (5)

This method has the advantage of being theoretically able to identify Pareto
optimal points also of non-convex problems. However, it also has two potential
drawbacks: the identified optimal point is only granted to be weakly Pareto
optimal, and the problem might become unfeasible due to the additional
constraints.

###  2.2 Goal Programming

The Goal Programming (or GP approach) is another widely used method to deal
with vector-valued problems (see [9] [10]). With respect to other
scalarization techniques, the idea behind this model is the determination of
the aspiration levels of an objective function. A GP model does not try to
find an optimal solution but an acceptable one, as it tries to achieve the
goals set by the DM rather than maximizing or minimizing the objective
functions. However, an optimization procedure is involved anyway. Within this
formulation, one tries to minimize any possible deviation from the objective
goals, either positive or negative. In fact, the GP model is a distance-
function model in which the obtained optimal solution represents the best
compromise between different objectives. Since the introduction of this
methodology for MOP problems, many variants have been presented in literature.
Among them, the most popular one is the Weighted Goal Programming ([1]) which
reads as follows: Given a set of ideal goals $g_{i}$, with $i=1,\dots,p$,
chosen by the DM, solve the following program:

$\min\sum_{i=1}^{p}\theta_{i}^{+}\delta_{i}^{+}+\theta_{i}^{-}\delta_{i}^{-}$

Subject to:

$\left\\{\begin{array}[]{l}J_{i}(x)+\delta_{i}^{-}-\delta_{i}^{+}=g_{i}\ \ \ \
i=1,\dots,p\\\ \delta_{i}^{-},\delta_{i}^{+}\geq 0\ \ \ \ \forall
i=1,\dots,p\\\ x\in X\end{array}\right.$ (6)

where $\delta_{i}^{+}$, $\delta_{i}^{-}$ are the positive and negative
deviations (slack variables), respectively, and $\theta_{i}^{+}$,
$\theta_{i}^{-}$ are the corresponding weights. Due to its simplicity, the GP
model and its variants have been widely applied to different areas such as
accounting, marketing, human resources, production, and so on. A negative
counterpart of the GP model, that it is important to consider, is the ability
of GP to produce solutions that are not Pareto efficient. To overcome this
difficulty, in order to produce Pareto optimal solutions the GP model is
implemented within a two-steps algorithm: At first, the GP solution is tested
for Pareto efficiency and, if it is not efficient, a restoration or projection
method is proposed to restore efficiency at the second step.

##  3 The Model

The simplest macroeconomic setting to account for intergenerational equity
([21]; [22]) and spatial links and spillovers ([6]; [3]) consists of an
optimization problem in which the social planner, by considering the evolution
of the capital stock $K(x,t)$, tries to determine the optimal level of
consumption $C(x,t)$ over time $t\in[0,T]$ and across space $x\in\Omega$ to
maximize the vector-valued social welfare which is composed of two terms. The
first term represents the discounted utility stream of the representative
individuals located in different venues within the entire spatial economy
$J_{1}=\int_{0}^{T}\int_{\Omega}U(C(x,t))e^{-\rho t}dxdt$, where the utility
function depends on consumption and the discount factor is $\rho>0$. The
second term represents the payoff of the representative individuals located in
different venues within the entire spatial economy at the end of the planning
horizon, $J_{2}=\int_{\Omega}K(x,T)$ where the payoff function depends on the
final capital level. The society’s optimization problem can therefore be
stated as a bi-criteria problem as follows:

$\max_{C(x,t)}W=\left(\int_{0}^{T}\int_{\Omega}U(C(x,t))e^{-\rho
t}dxdt,\int_{\Omega}K(x,T)dx\right)$ (7)

Subject to

$\small{\left\\{\begin{array}[]{ll}{\partial K(x,t)\over\partial
t}=\nabla(d(x)\nabla K(x,t))+f(K(x,t),C(x,t)),&(x,t)\in\Omega\times(0,T)\\\
d(x){\partial K\over\partial n}(x)=0,&x\in\partial\Omega\\\
K(x,0)=K_{0}(x).&x\in\Omega\end{array}\right.}$ (8)

where the spatio-temporal evolution of capital in (15) is driven by the
functions $f(K(x,t),C(x,t))$ with $f_{K}>0$ and $f_{C}<0$, and
$\nabla(d(x)\nabla K(x,t))$ which quantifies the process of spatial diffusion
between different locations within the spatial domain. Note that the function
$f(K(x,t),C(x,t))$ determines the nature of the trade-off between the two
terms in (7). The first term represents the traditional criterion considered
in macroeconomic analysis, characterizing the short termism of policymaking
which leads the society to care about the wellbeing of the current generation.
The second term represents instead a sustainability criterion introduced to
account for the wellbeing of future generations by allowing the society to
determine today how many resources not to exploit in order to preserve them
for the future. Clearly, a higher consumption today increases the first
criterion and through its subtraction of resources to capital accumulation
decreases the second one. Therefore, the society needs to balance these two
conflicting criteria determining the optimal dynamic path of consumption.

The above maximization has to be understood in a Pareto sense and with respect
to the Pareto order in ${\mathbb{R}}^{2}$: Given $a,b\in{\mathbb{R}}^{2}$,
$a\leq b$ if and only if $a_{i}\leq b_{i}$, $i=1,2$. So, in other words, a
feasible pair $(\bar{K}(x,t),\bar{C}(x,t))$ that solves (15) is optimal if
there is no other feasible pair $(K(x,t),C(x,t))$ such
that$(J_{1}(K(x,t),C(x,t)),J_{2}(K(x,t),C(x,t)))$ dominates
$(J_{1}(\bar{K}(x,t),\bar{C}(x,t)),J_{2}(\bar{K}(x,t),\bar{C}(x,t)))$. This
can be rewritten by stating that there exists no feasible $(K(x,t),C(x,t))$
such that

$(J_{1}(\bar{K},\bar{C}),J_{2}(\bar{K},\bar{C}))\in(J_{1}(K,C),J_{2}(K,C))\mathcal{}+{\rm
int}({\mathbb{R}}^{2}).$

Using the three approaches presented in the previous section to reduce a
multiple objective problem to a single objective model, we can define the
following three different single-criterion formulations.

###  3.1 Model I

In this first fomulation the two criteria are combined together through
scalarization weights. The DM decides the value of a trade-off parameter
$\Theta\in(0,1)$, that expresses the relative importance of first criterion
with respect to the second one. The parameter $\Theta$ is then expressing the
level of priority and preference that each criterion has for the DM. The
scalarized problem reads as

$\max J_{\Phi}(C,K):=(1-\Phi)J_{1}(C,K)+\Phi J_{2}(C,K)$

Subject to

$\small{\left\\{\begin{array}[]{ll}{\partial K(x,t)\over\partial
t}=\nabla(d(x)\nabla K(x,t))+f(K(x,t),C(x,t)),&(x,t)\in\Omega\times(0,T)\\\
d(x){\partial K\over\partial n}(x)=0,&x\in\partial\Omega\\\
K(x,0)=K_{0}(x).&x\in\Omega\end{array}\right.}$

where

$J_{1}(C,K)=\int_{0}^{T}\int_{\Omega}U(C(x,t))e^{-\rho t}dxdt$

and

$J_{2}(C,K)=\int_{\Omega}K(x,T)dx$

If we replace $\Theta={\Phi\over 1-\Phi}$ then the objective function boils
down to:

$\max J_{\Theta}(C,K):=J_{1}(C,K)+\Theta J_{2}(C,K)$

We will see in the following section dedicated to a numerical simulation that,
by varying over $\Theta\in(0,+\infty)$, it is possible to reconstruct the
Pareto frontier. In fact the objective function is concave in $(C,K)$ for any
$\Theta\in(0,+\infty)$.

###  3.2 Model II

In this second formulation we proceed by using the $\epsilon$-constraint
method. Within this approach there are essentially two different formulations
that can be proposed, once again this decision being dependent on the relative
importance of each criterion for the DM. If DM considers the intertemporal
utility as the main criterion and supposes that a certain level $\epsilon$ of
the criterion has to be attained, the model boils down to:

$\max J_{\epsilon}(C,K):=\int_{0}^{T}\int_{\Omega}U(C(x,t))e^{-\rho t}dxdt\\\
$ (9)

Subject to

$\small{\left\\{\begin{array}[]{ll}\int_{\Omega}K(x,T)dx\geq\epsilon,&\\\
{\partial K(x,t)\over\partial t}=\nabla(d(x)\nabla
K(x,t))+f(K(x,t),C(x,t)),&(x,t)\in\Omega\times(0,T)\\\ d(x){\partial
K\over\partial n}(x)=0,&x\in\partial\Omega\\\
K(x,0)=K_{0}(x).&x\in\Omega\end{array}\right.}$ (10)

By varying $\epsilon$ it is possible to obtain different points of the Pareto
frontier. If, instead, the DM is more interested in the sustainability
criterion, the model can be written as

$\max J_{\epsilon}(C,K):=\int_{\Omega}K(x,T)dx\\\ $ (11)

Subject to

$\small{\left\\{\begin{array}[]{ll}\int_{0}^{T}\int_{\Omega}U(C(x,t))e^{-\rho
t}dxdt\geq\epsilon,&\\\ {\partial K(x,t)\over\partial t}=\nabla(d(x)\nabla
K(x,t))+f(K(x,t),C(x,t)),&(x,t)\in\Omega\times(0,T)\\\ d(x){\partial
K\over\partial n}(x)=0,&x\in\partial\Omega\\\
K(x,0)=K_{0}(x).&x\in\Omega\end{array}\right.}$ (12)

###  3.3 Model III

In this third formulation we rely on a GP approach. Let us suppose that
$g_{1}$ and $g_{2}$ are the two goals associated with the two criteria $J_{1}$
and $J_{2}$, respectively. Within this technique, the DM is not interested in
the maximization process tout court, but instead in the achievement levels of
the two criteria $J_{1}$ and $J_{2}$. The previous model, reformulated using
the GP technique, reads as

$\min
J_{g_{1},g_{2}}:=\sum_{i=1}^{2}\theta_{i}^{+}\delta_{i}^{+}+\theta_{i}^{-}\delta_{i}^{-}$

Subject to

$\small{\left\\{\begin{array}[]{ll}\int_{\Omega}K(x,T)dx-\delta_{1}^{+}+\delta_{1}^{-}=g_{1}\\\
\int_{0}^{T}\int_{\Omega}U(C(x,t))e^{-\rho
t}dxdt-\delta_{2}^{+}+\delta_{2}^{-}=g_{2},&\\\ {\partial K(x,t)\over\partial
t}=\nabla(d(x)\nabla K(x,t))+f(K(x,t),C(x,t)),&(x,t)\in\Omega\times(0,T)\\\
d(x){\partial K\over\partial n}(x)=0,&x\in\partial\Omega\\\
K(x,0)=K_{0}(x).&x\in\Omega\end{array}\right.}$ (13)

##  4 Computational Results

In this section we present a numerical implementation of Model I introduced
above and we determine an approximation of the Pareto frontier. The next
section, instead, will present other possible numerical implementations of
Models II and III. Before doing that, we specify the components of this model
and, in particular, the form of the function $f(K(x,t),C(x,t))$. It is a
traditional Ramsey-type optimal control problem, extendeto account for
intergenerational equity issues ([12]; [13]) along with spatial heterogeneity
and spillovers ([3]; [16]). The economy develops along a linear city and the
social planner wishes to maximize the vector-valued social welfare by choosing
the level of consumption in each location, which in turn determines the
evolution of capital in each location and in the whole economy. Capital
accumulation depends on the difference between net (of depreciation, where the
depreciation rate is $\delta_{K}>0$) production, $Y(x,t)=AK(x,t)^{\alpha}$ and
consumption, augmented for the inflows of capital from other locations; these
flows are captured by a diffusion term $\nabla\left(d(x)\nabla K(x,t)\right)$,
where $d(x)$ is the diffusion parameter. Given the initial condition,
$K(x,0)=K_{0}(x)$, the problem can be summarized as follows:

$\displaystyle\max$
$\displaystyle\left(\int_{0}^{T}\int_{\Omega}U(C(x,t))e^{-\rho
t}dxdt,\int_{\Omega}K(x,T)dx\right)$ (14)

Subject to

$\small{\left\\{\begin{array}[]{ll}{\partial K(x,t)\over\partial
t}=\nabla\left(d(x)\nabla
K(x,t)\right)+AK(x,t)^{\alpha}-\delta_{K}K(x,t)-C(x,t),&(x,t)\in\Omega\times(0,T)\\\
d(x){\partial K\over\partial n}(x)=0,&x\in\partial\Omega\\\
K(x,0)=K_{0}(x).&x\in\Omega\end{array}\right.}$ (15)

We first discuss how to solve numerically the scalarized problem:

$\max J_{\Theta}(C,K):=J_{1}(C,K)+\Theta J_{2}(C,K)$

Subject to

$\small{\left\\{\begin{array}[]{ll}{\partial K(x,t)\over\partial
t}=\nabla\left(d(x)\nabla
K(x,t)\right)+AK(y,t)^{\alpha}dy-C(x,t)-\delta_{K}K(x,t),&(x,t)\in\Omega\times(0,T)\\\
d(x){\partial K\over\partial n}(x)=0,&x\in\partial\Omega\\\
K(x,0)=K_{0}(x).&x\in\Omega\end{array}\right.}$

To determine an optimal policy result, let us define the current Hamiltonian
function as

$H(C,K,\lambda)=U(C)+\lambda\left(\nabla(d(x)\nabla
K(x,t)+AK(x,t)^{\alpha}-\delta K(x,t)-C(x,t)\right)$

The following proposition provides the optimality conditions for an optimal
solution of the problem above.

######  Proposition 2.

Suppose that $U(C)$ is a concave function of $C$. Then a pair
$(\tilde{C},\tilde{K})$ solves the above optimal control model if and only if
it is solution to the following Hamiltonian system:

${\left\\{\begin{array}[]{ll}{\partial K(x,t)\over\partial
t}=\nabla\left(d(x)\nabla
K(x,t)\right)+AK(x,y)^{\alpha}-\delta_{K}K(x,t)-C(x,t),&(x,t)\in\Omega\times(0,T)\\\
{\partial\lambda(x,t)\over\partial
t}=\rho\lambda-\nabla\left(d(x)\nabla\lambda(x,t)\right)-\lambda\alpha
AK^{\alpha-1}(x,t)-\delta_{K}\lambda,&(x,t)\in\Omega\times(0,T)\\\
U^{\prime}(C)=\lambda&(x,t)\in\Omega\times(0,T)\\\ d(x){\partial
K\over\partial n}(x)=0,&x\in\partial\Omega\\\
d(x){\partial\lambda\over\partial n}(x)=0,&x\in\partial\Omega\\\
\lambda(x,T)=\Theta&x\in\Omega\\\
K(x,0)=K_{0}(x)&x\in\Omega\end{array}\right.}$

Since analyzing explicitly the Hamiltonian system above is generally not
possible (unless we introduce restrictive assumptions), we now proceed with
numerical simulations to illustrate the optimal behavior of capital and
consumption. If we use the dynamic constraint and plug it into the objective
function we obtain:

$J(K)=\int_{0}^{T}\int_{\Omega}U\left(-{\partial K(x,t)\over\partial
t}+\nabla\left(d_{k}(x)\nabla
K(x,t)\right)+AK(y,t)^{\alpha}-\delta_{K}K(x,t)\right)e^{-\rho
t}dxdt+\Theta\int_{\Omega}K(x,T)dx$

Subject to

$\left\\{\begin{array}[]{ll}-{\partial K(x,t)\over\partial
t}+\nabla\left(d(x)\nabla K(x,t)\right)+AK(x,t)^{\alpha}-\delta_{K}K(x,t)\geq
0,&(x,t)\in\Omega\times(0,T)\\\ d(x){\partial K\over\partial
n}(x)=0,&x\in\partial\Omega\\\ K(x,0)=K_{0}(x).&x\in\Omega\end{array}\right.$

The direction derivative of $J$ along any feasible $h$ is given by

$\displaystyle J^{\prime}(K;h)$ $\displaystyle=$ $\displaystyle\lim_{\delta\to
0}\frac{J(K+\delta h)-J(K)}{\delta}$ $\displaystyle=$
$\displaystyle\int_{0}^{T}\int_{\Omega}U^{\prime}(-{\partial
K(x,t)\over\partial t}+\nabla\left(d(x)\nabla
K(x,t)\right)+AK^{\alpha}(x,t)-\delta_{K}K(x,t))$ $\displaystyle*$
$\displaystyle\left(-{\partial h(x,t)\over\partial t}+\nabla\left(d(x)\nabla
h(x,t)\right)+Ah^{\alpha}(x,t)dx-\delta_{K}h(x,t)\right)e^{-\rho
t}dxdt+\Theta\int_{\Omega}h(x,T)dx$

We propose an algorithm to determine an approximation of the optimal solution.
At each step this algorithm determines the direction of growth $h$ using the
above calculated directional derivative $J^{\prime}(K;h)$.

  * •

Given the value of the state variable $K_{n}(x,t)$, solve the following
problem

${\left\\{\begin{array}[]{ll}-{\partial h(x,t)\over\partial
t}+\nabla\left(d(x)\nabla h(x,t)\right)+Ah^{\alpha}(x,t)-\delta_{K}h(x,t)=\\\
\left[U^{\prime}\left(-{\partial K_{n}(x,t)\over\partial
t}+\nabla\left(d(x)\nabla
K_{n}(x,t)\right)+AK_{n}^{\alpha}(x,t)-\delta_{K}K_{n}(x,t)\right)\right]^{-1}\left(-\Theta{\partial
h\over\partial t}+1\right)e^{\rho t},&(x,t)\in\Omega\times(0,T)\\\
d(x){\partial K\over\partial n}(x)=0,&x\in\partial\Omega\\\
h(x,0)=0.&x\in\Omega\end{array}\right.}$

  * •

Determine $\delta>0$ that corresponds to the maximum increment of $J$ along
the direction $h$

  * •

Update $K_{n+1}=K_{n}+\delta h$

  * •

If $|J(K_{n+1})-J(K_{n})|<\epsilon$ then stop otherwise go to the top.

The following result shows that $J$ is increasing along the sequence generated
by the above algorithm. The implementation of the above algorithm generates a
sequence of functions $K_{n}$ along which the objective function is
increasing.

######  Proposition 3.

If $\delta$ is small then $J(K_{n+1})\geq J(K_{n})$, $\forall n\geq 0$.

Proof. Computing we have:

$\displaystyle J(K_{n+1})-J(K_{n})$ $\displaystyle=$ $\displaystyle\delta
J^{\prime}(K_{n};h)+o(\delta)$ $\displaystyle=$
$\displaystyle\delta\left(\int_{0}^{T}\int_{\Omega}\left(-\Theta{\partial
h\over\partial t}+1\right)dxdt+\Theta\int_{\Omega}h(x,T)dx\right)+o(\delta)$
$\displaystyle=$
$\displaystyle\delta\left(-\Theta\int_{\Omega}h(x,T)-h(x,0)dx+T\mu(\Omega)+\Theta\int_{\Omega}h(x,T)dx\right)+o(\delta)$
$\displaystyle=$
$\displaystyle\delta\left(T\mu(\Omega)+{o(\delta)\over\delta}\right)\geq 0$

and this last passage relies on the boundary condition $h(x,0)=0$. ∎

We now apply the above algorithm to the following model where
$U(C)=[(1+C(x,t))^{2\over 3}-1]$, $T=1$, $\Omega=[0,1]$,
$d_{k}(x)=1-0.5x^{2}$, $\alpha=1$, $\delta_{K}=0.01$.

$\max\left(\int_{0}^{1}\int_{0}^{1}[(1+C(x,t))^{2\over
3}-1]dxdt,\int_{0}^{1}K(x,1)dx\right)$

Subject to

${\left\\{\begin{array}[]{ll}{\partial K(x,t)\over\partial
t}=\nabla\left((1-0.5x^{2})\nabla
K(x,t)\right)+K(x,t)-C(x,t)-0.01K(x,t),&(x,t)\in[0,1]\times(0,1)\\\
(1-0.5x^{2}){\partial K\over\partial n}(x)=0,&x=0,1\\\
K(x,0)=1+x.&x\in[0,1]\end{array}\right.}$

We have implemented the above algorithm in COMSOL. Firstly, we run the
algorithm with $\Theta=0$. The values of the objective function after three
iterations are: $J_{0}=0.579453721074241$, $J_{1}=0.6059822543917376$,
$J_{2}=0.6287663921318654$, $J_{3}=0.6534865860743782$. Then we used
$\Theta=0.1$. The values of the objective functions after three iterations are
given by: $J_{0}=0.5944537210742411$, $J_{1}=0.6200510708312351$,
$J_{2}=0.6404132122774402$, and $J_{3}=0.6570096551286734$. Figure 1 shows the
spatio-temporal behavior of $K$ in the $\Theta=0$ case (left panel) and
$\Theta=0.1$ case (right panel). It clearly shows our simple algorithm is
capable to handle an optimal control problem on partial differential equations
without requiring specific restrictive assumptions on the functional forms of
the utility and production functions.

Figure 1: Spatio-temporal evolution of capital in the $\Theta=0$ case (left)
and $\Theta=0.1$ case (right).

In order to construct the Pareto frontier we proceed as follows:

  * •

For each $\Theta$ we solve the “scalarized” problem $\max J_{1}+\Theta J_{2}$,

  * •

Let $K_{\Theta}(x,t)$ be the optimal solution, we plug it into the two
separated criteria and get the pair of values $(J_{1},J_{2})$,

  * •

We plot the pair $(J_{1},J_{2})$.

The Pareto frontier is illustrated in Figure 2, from which we can clearly
observe that it is bowed outward as expected, give the nature of the trade off
between the two criteria.

Figure 2: Piecewise linear Pareto frontier under scalarization.

##  5 Further Numerical Experiments

In this section we discuss some further numerical experiments that can be
implemented using different model formulations based on the
$\epsilon$-constraint and the GP approaches.

###  5.1 Experiment I

We now discuss a model formulation that has been obtained by applying the
$\epsilon$-constraint method. In this formulation we suppose that the DM
maximizes his intertemporal utility and the level of physical capital at the
final horizon $T$ is included in the set of constraints. Given a positive
value of $\epsilon$, the model can be written as:

$\max\int_{0}^{T}\int_{\Omega}U(C(x,t))e^{-\rho t}dxdt\\\ $ (16)

Subject to:

${\left\\{\begin{array}[]{ll}\int_{\Omega}K(x,T)dx\geq\epsilon,&\\\ {\partial
K(x,t)\over\partial t}=\nabla\left(d(x)\nabla
K(x,t)\right)+AK(x,t)^{\alpha}-\delta_{K}K(x,t)-C(x,t),&(x,t)\in\Omega\times(0,T)\\\
d(x){\partial K\over\partial n}(x)=0,&x\in\partial\Omega\\\
K(x,0)=K_{0}(x).&x\in\Omega\end{array}\right.}$ (17)

This model can be solved using the same approach and a slightly modified
version of the algorithm used in the linear scalarization case. If we use the
dynamic constraint to express the level of consumption $C(x,t)$ as function of
$K(x,t)$, the model boils down to:

$\max\int_{0}^{T}\int_{\Omega}U\left({\partial K(x,t)\over\partial
t}-\nabla\left(d(x)\nabla
K(x,t)\right)-AK(x,t)^{\alpha}+\delta_{K}K(x,t)\right)e^{-\rho t}dxdt\\\ $
(18)

Subject to:

${\left\\{\begin{array}[]{ll}\int_{\Omega}K(x,T)dx\geq\epsilon,&\\\ -{\partial
K(x,t)\over\partial t}+\nabla\left(d(x)\nabla
K(x,t)\right)+AK(x,t)^{\alpha}-\delta_{K}K(x,t)\geq
0,&(x,t)\in\Omega\times(0,T)\\\ d(x){\partial K\over\partial
n}(x)=0,&x\in\partial\Omega\\\ K(x,0)=K_{0}(x).&x\in\Omega\end{array}\right.}$
(19)

As a numerical experiment, we now apply the above model when
$U(C)=[(1+C(x,t))^{2\over 3}-1]$, $T=1$, $\Omega=[0,1]$,
$d_{k}(x)=1-0.5x^{2}$, $\alpha=1$, $\delta_{K}=0.01$. The model is formulated
as:

$\max\int_{0}^{1}\int_{0}^{1}\left[\left(1+-{\partial K(x,t)\over\partial
t}+\nabla\left((1-0.5x^{2})\nabla
K(x,t)\right)+K(x,t)-0.01K(x,t)\right)^{2\over 3}-1\right]dxdt$

Subject to

$\small{\left\\{\begin{array}[]{ll}\int_{0}^{1}K(x,1)dx\geq 1.3&\\\ -{\partial
K(x,t)\over\partial t}+\nabla\left((1-0.5x^{2})\nabla
K(x,t)\right)+K(x,t)-0.01K(x,t)\geq 0,&(x,t)\in[0,1]\times(0,1)\\\
(1-0.5x^{2}){\partial K\over\partial n}(x)=0,&x=0,1\\\
K(x,0)=1+x.&x\in[0,1]\end{array}\right.}$

The application of the above algorithm provides $J(1)=0.57945$,
$J(2)=0.60822$, $J(3)=0.63879$. The optimal behaviour of $K$ is shown in
Figure 3.

Figure 3: The evolution of $K$

###  5.2 Experiment II

In this model formulation we still use the $\epsilon$-constraint method but we
assume a linear utility function $U(C)=C$, and $\alpha=1$. In this context the
spatial model can be reduced to a one-dimensional optimal control model by
introducing the average of consumption and physical capital at $t$, $C_{M}(t)$
and $K_{M}(t)$, defined as

$C_{M}(t)=\int_{\Omega}C(x,t)dx$

and

$K_{M}(t)=\int_{\Omega}K(x,t)dx$

and the average amount of consumption per unit amount $c_{M}(t)$ as

$c_{M}(t)={C_{M}(t)\over K_{M}(t)}\in[0,1]$

Let us suppose that a pair $(K,C)$ solves the model:

$\max_{C(x,t),K(x,t)}\int_{0}^{T}\int_{\Omega}C(x,t)e^{-\rho t}dxdt\\\ $ (20)

Subject to:

${\left\\{\begin{array}[]{ll}\int_{\Omega}K(x,T)dx\geq\epsilon,&\\\ {\partial
K(x,t)\over\partial t}=\nabla\left(d(x)\nabla
K(x,t)\right)+AK(x,t)^{\alpha}-\delta_{K}K(x,t)-C(x,t),&(x,t)\in\Omega\times(0,T)\\\
d(x){\partial K\over\partial n}(x)=0,&x\in\partial\Omega\\\
K(x,0)=K_{0}(x).&x\in\Omega\end{array}\right.}$ (21)

Then, by easy calculations and taking the integral of the contraints, the pair
$(K_{M},c_{M})$ solves the following one:

$\max_{c_{M}(t),K_{M}(t)}\int_{0}^{T}c_{M}(t)K_{M}(t)e^{-\rho t}dt\\\ $ (22)

Subject to

${\left\\{\begin{array}[]{ll}\dot{K_{M}}(t)=AK_{M}(t)-\delta_{K}K_{M}(t)-c_{M}(t)K_{M}(t),&t\in(0,T)\\\
K_{M}(T)\geq\epsilon,&\\\
K_{M}(0)=\int_{\Omega}K_{0}(x)dx.&\end{array}\right.}$ (23)

This is classical calculus of variations model with a bang-bang control and it
can be solved in closed-form to determine the optimal paths of $C_{M}(t)$ and
$K_{M}(t)$.

###  5.3 Experiment III

This third experiment is devoted to a different formulation that has been
obtained by applying the GP approach. Here we also assume a linear utility
function $U(C)=C$, and $\alpha=1$. We also suppose that all weights are equal
and then normalized to $1$ and that $T$ is an integer number. Given two goals
$g_{1}$ and $g_{2}$ for the two criteria $J_{1}$ and $J_{2}$, the model reads
as:

$\min\delta_{1}^{+}+\delta_{1}^{-}+\delta_{2}^{+}+\delta_{2}^{-}$ (24)

Subject to:

${\left\\{\begin{array}[]{ll}\int_{\Omega}K(x,T)dx-\delta_{1}^{+}+\delta_{1}^{-}=g_{1}\\\
\int_{0}^{T}\int_{\Omega}C(x,t)e^{-\rho
t}dxdt-\delta_{2}^{+}+\delta_{2}^{-}=g_{2},&\\\ {\partial K(x,t)\over\partial
t}=\nabla\left(d(x)\nabla
K(x,t)\right)+AK(x,t)-\delta_{K}K(x,t)-C(x,t),&(x,t)\in\Omega\times(0,T)\\\
d(x){\partial K\over\partial n}(x)=0,&x\in\partial\Omega\\\
K(x,0)=K_{0}(x).&x\in\Omega\end{array}\right.}$ (25)

To simplify the model, we take the integral of the PDE and reduce the analysis
to the variables $K_{M}$ and $C_{M}$ introduced in the previous experiment.
The model boils down to:

$\min\delta_{1}^{+}+\delta_{1}^{-}+\delta_{2}^{+}+\delta_{2}^{-}$ (26)

Subject to:

${\left\\{\begin{array}[]{ll}K_{M}(T)-\delta_{1}^{+}+\delta_{1}^{-}=g_{1}\\\
\int_{0}^{T}C_{M}(t)e^{-\rho t}dt-\delta_{2}^{+}+\delta_{2}^{-}=g_{2},&\\\
\dot{K_{M}}(t)=(A-\delta_{K})K_{M}(t)-C_{M}(t),&t\in(0,T)\\\
K_{M}(0)=\int_{\Omega}K_{0}(x)dx.&\end{array}\right.}$ (27)

By discretizing the time, and introducing the discrete variables $K_{M}(i)$,
$C_{M}(i)$, $i=0...T$, the model can be written as

$\min\delta_{1}^{+}+\delta_{1}^{-}+\delta_{2}^{+}+\delta_{2}^{-}$ (28)

Subject to:

${\left\\{\begin{array}[]{ll}K_{M}(T)-\delta_{1}^{+}+\delta_{1}^{-}=g_{1}\\\
\sum_{i=0}^{T}C_{M}(i)e^{-\rho i}-\delta_{2}^{+}+\delta_{2}^{-}=g_{2},&\\\
K_{M}(i+1)-K_{M}(i)=(A-\delta_{K})K_{M}(i)-C_{M}(i),&i=0...T-1\\\
K_{M}(0)=\int_{\Omega}K_{0}(x)dx.&\end{array}\right.}$ (29)

This is a linear optimization model that can be solved by standard
optimization solvers such LINGO or MATLAB.

##  6 Conclusion

The introduction of intergenerational equity considerations associated with
sustainability issues into a traditional macroeconomic setting transforms the
typical unicriterion macroeconomic problem into a bi-criteria optimization
problem, which can be analyzed through the lenses of the multicriteria
optimization techniques developed in the operations research literature.
Recently traditional macroeconomic problems have been extended to introduce a
spatial dimension, allowing to consider the extent to which spatial
heterogeneity and spatial spillovers affect economic outcomes not only over
time but also across space. The goal of this papers consists thus to merge
these two different lines of research by analyzing a simple macroeconomic
setting to account for intergenerational equity and spatial spillovers from
the operations research point of view. In particular, we show that our
macroeconomic problem can be reformulated as a multicriteria problem by
relying on different techniques (scalarization, $\epsilon$-constraint method
and goal programming), and such different formulations of the problem can be
solved through numerical methods, which allow us to illustrate the nature of
the trade-off between the two criteria and to derive the Pareto-frontier in
some specific cases and parametrizations.

Our paper represents one of the first attempts to bridge the economics and the
operations research literature, but still much needs to be done in order to
develop further the possible synergies existing between these two different
disciplines. In particular, apart from scalarization, $\epsilon$-constraint
method and goal programming approaches, several other techniques developed in
the operations research literature can be applied in similar macroeconomic
contexts, especially in the context of stochastic or fuzzy multiple objective
optimization. It would also be worth exploring the use of vectorization
algorithms and methods that tackle the MOP model directly such as, for
instance, genetic algorithms. Moreover, apart from the applications to
macroeconomic questions, similar multicriteria approaches can be applied to
other economic problems arising in environmental economics, game theory and
cost-benefit analysis.

## References

  * [1] Aouni B., Colapinto C. and La Torre (2014). Financial portfolio management through the goal programming model: Current state-of-the-art, European Journal of Operational Research 234, pp. 536-545.
  * [2] Arrow, K., Dasgupta, P., Goulder, L., Daily, G., Ehrlich, P., Heal, G., Levin, S., Maler, K.G., Schneider, S., Starrett D., Walker, B. (2004). Are we consuming too much?, Journal of Economic Perspectives 18, 147–172
  * [3] Boucekkine, R., Camacho, C., Zou, B. (2009). Bridging the gap between growth theory and economic geography: the spatial Ramsey model, Macroeconomic Dynamics 13, 20–45
  * [4] Boucekkine, R., Camacho, C., Fabbri, G. (2013a). On the optimal control of some parabolic differential equations arising in economics, Serdica Mathematical Journal 39, 331–354
  * [5] Boucekkine, R., Camacho, C., Fabbri, G. (2013b). Spatial dynamics and convergence: the spatial AK model, Journal of Economic Theory 148, 2719–2736
  * [6] Brito, P. (2004). The dynamics of growth and distribution in a spatially heterogeneous world, UECE-ISEG, Technical University of Lisbon
  * [7] Camacho, C., Zou, B. (2004). The spatial Solow model, Economics Bulletin 18, 1–11
  * [8] Camacho, C., Zou, B., Briani, M. (2008). On the dynamics of capital accumulation across space, European Journal of Operational Research 186 2, 451–465
  * [9] Charnes, A., Cooper, W.W., Ferguson, R.O. (1955). Optimal estimation of executive compensation by linear programming, Management Science, 2 138–151
  * [10] Charnes, A., Cooper, W.W. (1961). Management models and industrial applications of linear programming (New York: John Wiley & Sons)
  * [11] Chinchilnisky, G., Heal, G., Beltratti, A. (1995). The green golden rule, Economics Letters 49, 174–179
  * [12] Chinchilnisky, G. (1997). What is sustainable development?, Land Economics 73, 476–491
  * [13] Colapinto, C., Liuzzi, D., Marsiglio, S. (2017). Sustainability and intertemporal equity: a multicriteria approach, Annals of Operations Research 251, 271–284
  * [14] Heal, G. (2005). Intertemporal welfare economics and the environment, in (Maler, K.G., Vincent, J.R., Eds.), “Handbook of Environmental Economics”, vol. 3 (North-Holland: Amsterdam)
  * [15] Krugman, P. (1991). Increasing returns and economic geography, Journal of Political Economy 99, 483–499
  * [16] La Torre, D., Liuzzi, D., Marsiglio, S. (2015). Pollution diffusion and abatement activities across space and over time, Mathematical Social Sciences 78, 48–63
  * [17] La Torre, D., Liuzzi, D., Marsiglio, S. (2017). Pollution control under uncertainty and sustainability concern, Environmental and Resource Economics 67, 885–-903
  * [18] La Torre, D., Liuzzi, D., Marsiglio, S. (2019). Population and geography do matter for sustainable development, Environment and Development Economics 24, 201-–223
  * [19] La Torre, D., Marsiglio, S. (2019). A note on optimal debt reduction policies, Macroeconomic Dynamics, forthcoming
  * [20] Le Kama, A.D.A. (2001). Sustainable growth, renewable resources and pollution, Journal of Economic Dynamics & Control 25, 1911–1918
  * [21] Marsiglio, S., La Torre, D. (2018). Economic growth and abatement activities in a stochastic environment: a multi-objective approach, Annals of Operations Research 267, 321-–334
  * [22] Marsiglio, S., Privileggi, F. (2019). On the economic growth and environmental trade–off: a multi-objective analysis, Annals of Operations Research, forthcoming
  * [23] Pezzey, J.C.V. (1997). Sustainability constraints versus “optimality” versus intertemporal concern, and axioms versus data, Land Economics 73, 448–466
  * [24] Ramsey, F. (1928). A mathematical theory of saving, Economic Journal 38, 543–559
  * [25] Quah, D.T. (1996). Regional convergence clusters across Europe, European Economic Review 40, 951–958
  * [26] Ramsey, F. (1928). A mathematical theory of saving, Economic Journal 38, 543–559
  * [27] Sawaragi, Y., Nakayama, H., Tanino, T. (1985). Theory of multiobjective optimization (Academic Press, Inc.)
  * [28] Solow, R.M. (1956). A contribution to the theory of economic growth, Quarterly Journal of Economics 70, 65–94
  * [29] Wolfe, D.A., Gertler, M.S. (2004). Clusters from the inside and out: local dynamics and global linkages, Urban Studies 41, 1071–-1093
  * [30] World Commission on Environment and Development (1987). Our common future (Oxford University Press, Oxford)