Mathematical modeling and optimal control of carbon dioxide emissions from energy sector

Abstract

Energy demand is rising day by day and will continue to increase to meet the demand of the growing population. A major portion of global energy production comes from fossil fuel burning, resulting in the increase in the atmospheric burden of global warming gas carbon dioxide ($CO _{2}$). Cutting down $CO _{2}$ emission from the energy sector is crucial to meet the climate change mitigation target. This paper is focused on fulfilling two objectives: The first objective is to present a mathematical model that captures the dynamical relationship between the human population, energy use, and atmospheric carbon dioxide, and the second aim is to derive a mathematical framework to effectively utilize the available mitigation options to curtail $CO _{2}$ emission from energy use by proposing an optimal control problem. The mitigation options that reduce the $CO _{2}$ emission rate from energy production, as well as the options that reduce the energy consumption rate, are considered in the modeling process. The proposed mathematical model is analyzed qualitatively to comprehend the system’s long-term behavior. The model parameters are fitted to real data of global energy use, population, and $CO _{2}$ concentration. It is shown that the equilibrium level of $CO _{2}$ reduces with the increase in the efficiencies of mitigation options to reduce the $CO _{2}$ emission rate per unit energy use and energy consumption rate. The optimality system is derived analytically by taking the efficiencies of the mitigation options to reduce the $CO _{2}$ emission rate and energy consumption rate as control variables. Numerical simulations are conducted to validate the theoretical findings and identify the optimal profiles of control variables under different settings of $CO _{2}$ emission rate, energy consumption rate, and maximum efficiencies of available mitigation options to cut down $CO _{2}$ emission rate and energy consumption rate. It is found that the development and implementation of more efficient mitigation options and switching to low carbon energy sources bring reduction in the mitigation cost.

Appendices

Appendix

A Proof of Theorem 1

Let $J_{S_i}$ denote the Jacobian matrix of system (4) evaluated at $S_i$. The eigenvalues of $J_{S_1}$ are $-\alpha$, r, and 0. As one eigenvalue is $r>0$, therefore $S_1$ is unstable.

One of the eigenvalues of $J_{S_2}$ is $\frac{\gamma (1-\mu _1) N_2}{K+N_2}$. The other two eigenvalues are root of equation $x^2+\left( \alpha +\frac{r N_2}{L}\right) x+ \left( \frac{\alpha r}{L}+\theta \lambda _1\right) N_2=0$, which are either negative or with negative real part. Since one eigenvalue of $J_{S_2}$ is positive, therefore $S_2$ is also unstable.

The local stability of $S^*$ is examined using Lyapunov’s direct method. Consider the following function:

$$\begin{aligned} W=\frac{1}{2}c^2+\frac{p_1}{2} \frac{n^2}{N^*}+\frac{p_2}{2} \frac{e^2}{E^*}, \end{aligned}$$

(23)

where $p_1$ and $p_2$ are positive constants. Here c, n and e are small perturbations in C, N and E about the steady state $S^*$.

The function W is positive definite function. The time derivative of ‘W’ along the linearized system of (4) corresponding to $S^*$ is given by

$$\begin{aligned} \dot{W}= & {} -\alpha c^2-p_1\left( \frac{r}{L}-\beta _2E^*\right) n^2-p_2\gamma _0 e^2+(\lambda _1-p_1 \theta )c n\\&+\lambda _2(1-\mu _2) c \ e+p_1(\beta _1 +\beta _2N^*)n \ e+\frac{p_2(1-\mu _1)\gamma K}{(K+N^*)^2}n \ e. \end{aligned}$$

After taking $p_1=\frac{\lambda _1}{\theta }$, $\dot{W}$ is negative definite if the following inequalities hold:

$$\begin{aligned}&p_2 >\frac{3}{4}\frac{\lambda _2^2(1-\mu _2)^2}{\alpha \gamma _0} \end{aligned}$$

(24)

$$\begin{aligned}&p_2 > \frac{3}{2}\frac{\lambda _1(\beta _1+\beta _2 N^*)^2}{\theta \left( \frac{r}{L}-\beta _2 E^*\right) \gamma _0} \end{aligned}$$

(25)

$$\begin{aligned}&p_2 < \frac{2}{3}\frac{\lambda _1 \gamma _0(K+N^*)^4\left( \frac{r}{L}-\beta _2 E^*\right) }{\theta (1-\mu _1)^2K^2\gamma ^2} \end{aligned}$$

(26)

From the above inequalities, it is found that $\dot{W}$ is negative definite, and hence, W is a Lyapunov function, under the condition (12).

B Proof of Theorem 2

Consider the following scalar valued positive definite function:

$$\begin{aligned} V=\frac{1}{2}(C-C^*)^2+m_1\left( N-N^*-N^* ln \frac{N}{N^*}\right) +m_2\left( E-E^*-E^*ln \frac{E}{E^*}\right) , \end{aligned}$$

(27)

where $m_1$ and $m_2$ are positive constants.

The time derivative of ‘V’ is given as

$$\begin{aligned} \dot{V}= & {} -\alpha (C-C^*)^2-m_1\left( \frac{r}{L}-\beta _2E^*\right) (N-N^*)^2-m_2\gamma _0(E-E^*)^2\\&+(\lambda _1-m_1 \theta )(C-C^*)(N-N^*)+\lambda _2(1-\mu _2)(C-C^*)(E-E^*)\\&+m_1(\beta _1 +\beta _2N)(N-N^*)(E-E^*)+\frac{m_2(1-\mu _1)\gamma K}{(K+N)(K+N^*)}(N-N^*)(E-E^*). \end{aligned}$$

Choosing $m_1=\frac{\lambda _1}{\theta }$, $\dot{V}$ is negative definite if the following inequalities hold:

$$\begin{aligned}&m_2 >\frac{3}{4}\frac{\lambda _2^2(1-\mu _2)^2}{\alpha \gamma _0} \end{aligned}$$

(28)

$$\begin{aligned}&m_2 > \frac{3}{2}\frac{\lambda _1(\beta _1+\beta _2 N_m)^2}{\theta \left( \frac{r}{L}-\beta _2 E^*\right) \gamma _0} \end{aligned}$$

(29)

$$\begin{aligned}&m_2 < \frac{2}{3}\frac{\lambda _1 \gamma _0(K+N^*)^2\left( \frac{r}{L}-\beta _2 E^*\right) }{\theta (1-\mu _1)^2\gamma ^2} \end{aligned}$$

(30)

The above inequalities are reduced to the condition (13). Thus, V is Lyapunov function on $\varOmega$ provided the condition (13) holds.