Global warming has been one of the most crucial problems, and maybe the greatest related to our real-world environment, which humankind is seriously faced with. As a result of industrialization and economic development, emission of greenhouse gases, specially carbon dioxide (CO 2 ), caused a significant increase in the mean temperature of the atmosphere by 1 • C to 1.5 • C, according to international reports [1], relative to pre-industrial era. This deviation of temperature, caused many changes in the global climate and environment, including but not limited to melting of ice lands [2], increasing trend towards larger wildfires [3], [4], increased heat content of deep level oceans [5], sea level rise [6], increased heavy precipitation [7], and earlier timing of spring events [2], [8]. According to a special report [9] from Inter-governmental Panel on Climate Change (IPCC), which is published in 2018, to keep the temperature deviation below 2 • C, we will require to reach zero greenhouse gases emission by year 2050. Despite the fact that this target temperature is remaining in very same condition of present state, however, reaching this objective is very hard and it needs a global cooperation among all countries and nations. In order to prevent the destructive consequences of climate change, it must be taken care of every little bit of emission. Scientists of the field use the Integrated Assessment Models (IAMs) to model the economy-climate interactions and quantify the damages of greenhouse emissions on the society and life on the planet. While many models proposed to describe the climate change economy, only three of them are mostly cited and commonly used as IAMs, all of which make possible estimation of Social Cost of Carbon (SCC). These models are as follows: Dynamic Integrated model of Climate and Economy (DICE) [10]- [17], Climate Framework for Uncertainty, Negotiation, and Distribution (FUND) [18], and Policy Analysis for the Greenhouse Effect (PAGE) [19]. As an example, the United States Interagency Working Group (IWC) used these three models to estimate the social cost of carbon dioxide emissions. These three models, DICE, FUND, and PAGE, are different. For example, while DICE is a general equilibrium model, as it models the economic growth as a state variable in the model, PAGE and FUND models are partial equilibrium models, hence, they do need to have an external input describing the economic factors. Also for FUND and PAGE models, there is not any optimization problems associated with the model, therefore, it is impossible to define an optimal control problem on these models. However, DICE is associated with an objective function and it can be expressed as an optimal control problem, accordingly. In this paper, we are going to use DICE and define a multiobjective optimal control problem on this model. Perhaps, Professor William Nordhaus proposed first IAM model of the field in [10] and then developed that model to reach the DICE. First version of DICE is proposed by Nordhaus in 1992 [11], and he improved and extended the model regularly in several years [12]- [17]. Alongside the DICE model, there is a regional version of this model, named Regional Integrated model of Climate and Economy (RICE) [14], which divides the whole world into 12 economic regions with their specific economic dynamics. In 2018, for his work in the field of climate change economy, Nordhaus received the Nobel Memorial Prize in Economic Sciences. There are two versions of DICE model, which are mostly cited in the literature: DICE2013 [15] and DICE2016 [16], [17]. These two models have not any structural differences. But some of parameters and time window of these two editions of DICE model are different. In this paper, we will use the dynamics and parameters of the DICE2016 version. Previously, Model Predictive Control (MPC) is used to solve the optimal control problem on DICE model in [20] and [21]. In these works, as well as original works of Nordhaus, the objective function is defined to be the Social Welfare, which incorporates both economic and environmental factors. However, to have a direct approach towards optimization of climate change and temperature deviation, while keeping optimizing the economic factors, we converted the problem into a multi-objective optimal control problem and used Non-Dominated Sorting Genetic Algorithm II (NSGA-II) to solve it. The results indicate that not only the result of MPC [20], [21] is discoverable by NSGA-II, but by having a complete set of Pareto-optimal solutions, we can make better decisions to reach lower temperature deviations levels. This paper is organized as follows. In Section II we will describe the dynamics and structure of DICE model in detail. In Section III, we will define a multi-objective optimal control problem on DICE model, which we are going to solve it, in following sections. Section IV contains a brief description of NSGA-II, as well as the structure of solutions to the problem and the method used for generating initial population. In Section V, the results of the NSGA-II for the problem are presented and finally, Section VI concludes the paper. II. DICE MODEL DICE model is a discrete-time model, and its latest version, starts from year 2015 as initial time. The time and discrete index are related to each other, according to the following equation:

where i, t, t 0 and ∆t denote to discrete time index, time, starting year (2015), and sample time (5 years), respectively. We will use the discrete index in the equations of DICE model. Total population of the world is given by the following dynamic equation:

where L(i) denotes the global population at time step i, and L a and g are parameters. Equation (2) indicates that dynamics of global population is assumed to be independent from any other state variables of the system. The parameter L a is the maximum value of worldwide population L(i). For current time, the exponential term in right-hand side of (2) is greater than 1, because L(i) L a for now. This means the population will increase by time.
However, as global population increases, the exponential term will be closer and closer to 1. This will keep the population below L a and just in limit, L(i) will reach it. According to Cobb-Douglas production function, the Gross Economic Output is expressed as:

where A(i), K(i) and γ denote technological efficiency (aka. total factor productivity), total capital, and capital elasticity constant, respectively. Technological efficiency A(i) is driven by the following dynamics:

where g A and δ A are constant parameters. This equation expresses the dynamics of another state variable A(i), which is assumed to be independent from other states in the DICE model. Denominator of the right-hand side of (4) will converge to 1 -g A as i → ∞. Therefore, in limit, the total progress of technological efficiency within a 5-year period, equals to 1 1-g A , which is approximately equal to 1.0823. This means, according to DICE model, technological efficiency improves about 8.23% within 5 years, at most. Total capital K(i)is given by the following dynamical equation:

where δ K and I(i) are respectively annual capital depreciation rate and total annual investment at time step i. Equation ( 5) says that, either capital is invested or remains as capital and depreciates. The amount of investment is defined as:

where s(i) ∈ [0, 1] is a control input, known as saving rate, and Q(i) is Net Economic Output. After subtracting the investment from net economic output, the remaining amount will be total consumption C(i), which is given by:

Net economic output is defined as follows:

where Λ(i) and Ω(i) are abatement cost fraction and climate damage function, respectively. These functions reduce the gross economic output Y (i) to get the net output Q(i) and they are defined by following relations:

and

where µ(i) ∈ [0, 1] is a control input, known as mitigation rate or abatement fraction, θ 1 (i) is the cost of mitigation effort, T AT (i) is deviation of atmosphere temperature with respect to year 1900, and θ 2 , ψ 1 , and ψ 2 are parameters. Dynamics for T AT (i) will be discussed in the next section. In damage function, ψ 1 is set to be zero and ψ 2 is set to have a total damage of 2 percent for T AT (i) = 3 • C. Also, the cost of mitigation effort θ 1 (i) is given by the following equation:

where σ(i) is the global emission intensity of economic activity, which is defined in the next section, and p b , δ pb , and θ 2 are parameters. The mitigation or abatement cost fraction Λ(i) is related to emission intensity, hence, it connects the net economic output to carbon emission dynamics. On the other hand, the damage function Ω(i), is related to T AT (i), i.e. atmosphere temperature deviation. This relates the economic dynamics to climate and carbon emissions in advance. Actually, ( 8) is the part of model which relates the socio-economic parts of the model to the environment and climate. Global emission intensity dynamics is governed by:

where g σ and δ σ are parameters. This equation indicates that emission intensity only depends on time and it is not related to other state variables of the mode. The equilibrium point of ( 12) is σ(∞) = 0. Because denominator of the right-hand side of the equation is always greater than 1. Hence, the value of emission intensity will be always decreasing, until it reaches zero in limit. Total emission is the sum of emissions due to economic activity and emissions related to the change of land use and it is given by:

where E(i) and E Land (i) are total emissions and land-userelated emissions, respectively. Equation ( 13) says that per one unit of gross economic output Y (i), total emissions increase by σ(i) units. Here, the mitigation control input µ(i) affects the emissions, and in case of full mitigation, total amount of economy related emissions can be cancelled out. Emissions related to the change of land use is given by:

where E L,0 is the estimation of initial emission due to land use in base year (2015), and δ EL is respective reduction rate. This equation indicates that land-use-related emission only depends on time and it is reducing as time passes. In this model, the total possible amount of emissions related to land use, from now to far future, is assumed to be constant, and according to (14) it can be calculated as (1 + 1 δEL )E L,0 . In order to model the mass of carbon on the planet, three state variables introduced: the average mass of carbon in atmosphere, upper ocean carbon mass, and lower (deep) ocean carbon mass, which are denoted by M AT , M UP , M LO , respectively. Dynamics of these quantities are given by:

where the elements of state transition and input matrices of this state-space model are constant parameters. According to (15), carbon mass in atmosphere is affected by previous values of atmosphere and upper ocean carbon mass, using constants parameters ζ 11 and ζ 12 , respectively. Similarly, upper ocean carbon mass is related the past values of carbon mass in atmosphere, upper ocean and lower ocean. Finally, deep ocean carbon mass is not affected by carbon mass in atmosphere and it is affected by past values of ocean related carbon mass, both upper and lower ocean. Also, total emissions E(i) only increases the amount of carbon mass in atmosphere, and according to this model, it takes more than 5 years to see the effect of emissions in carbon mass of upper or lower ocean levels. Temperature deviation of atmosphere, denoted by T AT (i), and lower (deep) ocean temperature deviation, denoted by T LO (i), are two states to describe the climate model in DICE. These are related to each other and they are driven by the following dynamics: (16) where elements of state transition and input matrices are constant parameters.
The state transition matrix is full, which means both atmosphere and deep ocean temperature deviations, i.e. T AT (i) and T LO (i), affect each other. The input to temperature dynamics, F (i), is radiative force resulted by greenhouse effect, which obviously should alter the atmosphere temperature and it is given by:

where F 2× is the force related to doubling of atmospheric carbon mass, M AT,1750 is the carbon mass in atmosphere in year 1750, and F EX (i) is an exogenous force, described as:

with parameters f 0 , f 1 and t f as minimum force, maximum force, and slope factor, respectively. The starting value of this exogenous force is f 0 , and it increases by 1 t f every time step, i.e. 5 years, until it reaches the f 1 . After that, it remains constant. Maximization of social welfare function is the original objective of DICE model. The objective function is defined as a discounted sum of utilities within an infinite horizon and it is given by:

where ρ is the discount rate, C(i) is total consumption, L(i) is total population, and utility function U (C, L) is defined as:

with α ≥ 0 as elasticity of marginal utility of consumption. Maximizing the value of W (µ, s), subject to dynamics of DICE model is a nonlinear optimal control problem. Because utility function given in (20) and social welfare defined by ( 19) are related to consumption, according to ( 7) and ( 8), the objective function is related to both socio-economic factors and climate conditions. However, the function reflects the economic part more than climate-related factors. Hence, in this paper, we are going to introduce another objective, which is directly related to climate. Parameter values of DICE model are gathered in Table I. These values are updated in [16], [17] and used in DICE-2016R version of the model. In this paper, we are going to modify the original optimal control problem on DICE model into a multi-objective problem. Actually the welfare function defined by ( 19) contains factors from socioeconomic and climate dynamics, and it reflects the optimality in both directions. However, having an additional straight objective related to climate or carbon mass will provide a complete set of options which can lead to better decisions. We define a bi-objective optimal control problem, to maximize social welfare and minimize the maximum atmosphere temperature deviation within a specified time window:

subject to dynamics provided by DICE model, i.e. dynamic equations ( 2), ( 4), ( 5), ( 12), ( 15) and ( 16), and assuming control inputs 0 ≤ µ(i) ≤ 1 and 0 ≤ s(i) ≤ 1, for every time index i. Time horizon of this problem is H; hence, i varies from 0 to H. Therefor, T AT,max is defined by following relation:

Our second objective is directly related to deviation in temperature and it is a kind of worst-case optimization criterion, very similar to what used in robust control, e.g. H ∞ control. Solving the multi-objective control problem, defined by ( 21), will results in various options, known as Pareto solutions, and this makes possible to make sacrifices in social welfare for the sake of lower temperature deviation and better climate. But single-objective optimal control, misses all of these possibilities and stocks in a specific solution, without any flexibilities. It is also possible to choose other objectives, for example maximizing the net economic output while minimizing the mass of atmosphere carbon. But according to our tests, the combination of welfare function and maximum temperature deviations, has much better results and provides diverse and meaningful options. In this paper, the multi-objective optimal control defined in previous section, is solved using Non-Dominated Sorting Genetic Algorithm (NSGA-II). In this section, we will briefly review the structure of this algorithm and the approached we used to solve the problem. Originally, Srinivas and Deb [22] proposed a new kind of genetic algorithm to deal with multi-objective optimization problems and called it Non-Dominated Sorting Genetic Algorithm (NSGA), which uses the multi-objective domination as ranking (sorting) criterion and fitness sharing mechanism as diversification mechanism. Performance of NSGA is highly sensitive to parameters of fitness sharing procedure. Hence, Deb et al. [23] introduced a new version of the algorithm, named NSGA-II, which is much faster and less sensitive to parameters. Instead of fitness sharing, NSGA-II uses crowding distance to control the diversity of solutions. After initialization, at every iteration of NSGA-II, the crossover and mutation operations are applied to selected members of current population to create a population of offsprings and mutants Usually, parents are selected using binary tournament selection mechanism in NSGA-II. After creation of new individuals, they are merged into current main population to create a larger population. The resulting population, is sorted according to two criteria: (a) non-domination rank, which is determined by non-dominated sorting procedure, and (b) crowding distance, which controls the diversification by removing similar solutions and keeping distinct ones. Finally, members of the population with rank 1 (first front, F 1 ), will form Pareto frontier. A pseudo-code of NSGA-II is given in Fig. 1. Solutions returned by NSGA-II are not dominated by other solutions, i.e. they are first-rank members of the population. These solutions form Pareto front and they can be seen as multi-objective optimal, or at least sub-optimal, solutions of

According to ( 21) and ( 22), we have two control inputs and we must determine the value of these signals at every time step. So the total number of decision variables is 2 × H, and all of them are in the range [0, 1]. The structure of the solutions is as follows:

(23) Note that input at time step H will not have any effects on the states, which are used to calculate the objective functions. Actually, they affect the states at time step H + 1, which is out of problem time window.
In the initialization phase of the algorithm, solutions are created assuming constant inputs over time, i.e. all elements in first row of solution matrix have same value, as well as the elements of second row, sharing their own common value. Creating initial solutions according to this scheme, results in a better Pareto frontier. and the time horizon (H) is set to be 37. So we seek for Paretooptimal values of decision variables in a 74-dimensional search space. The time window of our simulation of DICE model, is from year 2015 to 2200. Fig. 2 illustrates the resulting Pareto frontier. To compare the results with [20] and [21], the solution of Model Predictive Control (MPC) is shown in bi-objective plane, as well. According to this plot, result of MPC is very close to be Pareto-optimal. However, the solution provided by MPC is dominated by some points on the Pareto frontier and it is not a Pareto-optimal solution. Having this Pareto front, makes possible to make better decisions, for example compromise about reducing in socioeconomic objective of welfare, to have better climate. Some of solutions, marked with A to F, are selected from Pareto set, with equal spacing with respect to temperature deviation (first objective). Solution A is related to maximum social welfare (best) and maximum temperature deviation (worst). Similarly, solution F is related to minimum deviation of atmosphere temperature (best) and minimum welfare (worst). According to path from F to A, increasing social welfare will result the temperature deviation, and this is because of conflicting nature of objective functions. As mentioned earlier in subsection II-E, social welfare function is related to both socio-economic and climate-related factors. However, it is mainly affected by economic variables, such as gross economic output and it is less sensitive to environment conditions. This can be inferred from Fig. 2, too. The result of MPC is very close to the upper-left part of the Pareto frontier. The single-objective approach, which is aimed to maximize the social welfare only, discovers a solution with almost highest possible welfare and temperature deviation. This show the true nature of social welfare function and its sensitivity to climate-related factors, such as atmosphere temperature deviation. Some of states and control inputs for these selected solutions are shown in Fig. 3 to Fig. 6. These plots illustrate the deviation of atmosphere temperature (from temperature in year 1900) T AT , total emissions E, mitigation rate µ (first control input), and saving rate s (second control input), respectively. In addition to selected solutions from Pareto set, the same signals are plotted for MPC results [20], [21], too. Objective values of selected solutions, A to F, and solution found by MPC are listed in Table II. The rows of table are sorted in descending order of social welfare value. According to the objective values, for example, moving from MPC to D, will reduce the welfare by 0.0077 (less than 0.03%), but it will reduce the temperature deviation by more than 1.15 • C (more than 26 percent, relatively). However, the results indicate that it is not possible to have temperature deviations less than 2.37 • C, with current technologies and efforts. In best situation, this is the minimum reachable temperature, if everything goes as predicted by DICE model. Unless, some technological advancements will emerge, and makes it possible to go below this limit. A bi-objective optimal control problem is defined on DICE model, to find Pareto-optimal strategies for carbon mitigation (abatement) and saving rate (investment), which maximize the social welfare function while minimizing the deviation of atmosphere temperature. This problem is solved by NSGA-II, a multi-objective genetic algorithm and a Pareto frontier is found. There are solutions in this set which dominate the result of previous works on model predictive control of DICE model. Also, the Pareto front offers a wide range of options to reach lower temperature deviation by sacrificing welfare. However, the results show that, unfortunately, it is not possible to go below some limits and decrease the temperature deviation as much as desired. Application of the same technique on RICE model, a regional version of DICE which is also developed by Nordhaus, and changing the structure of algorithm to find smooth control inputs, can be possible directions for future research in this field.