text:
\begin{document}

\def\spacingset#1{\renewcommand{\baselinestretch}
{#1}\small\normalsize} \spacingset{1}

  \title{\bf Discovering an Active Subspace in a\\ Single-Diode Solar Cell Model}
  \author{Paul G.~Constantine\thanks{
    The authors gratefully acknowledge the support Department of Energy's Advanced Scientific Computing Research Award \emph{Mathematical and Statistical Methodologies for DOE Data-Centric Science at Scale}.}\hspace{.2cm}
    and
    Brian Zaharatos\\
    Department of Applied Mathematics and Statistics\\
    Colorado School of Mines\\
    and \\
    Mark Campanelli\\
    National Renewable Energy Laboratory}
  \maketitle

\bigskip
\begin{abstract}
Predictions from science and engineering models depend on the values of the model's input parameters. As the number of parameters increases, algorithmic parameter studies like optimization or uncertainty quantification require many more model evaluations. One way to combat this curse of dimensionality is to seek an alternative parameterization with fewer variables that produces comparable predictions. The \emph{active subspace} is a low-dimensional linear subspace defined by important directions in the model's input space; input perturbations along these directions change the model's prediction more, on average, than perturbations orthogonal to the important directions. We describe a method for checking if a model admits an exploitable active subspace, and we apply this method to a single-diode solar cell model with five input parameters. We find that the maximum power of the solar cell has a dominant one-dimensional active subspace, which enables us to perform thorough parameter studies in one dimension instead of five.
\end{abstract}

\noindent
{\it Keywords:} single-diode solar cell model, active subspaces, dimension reduction, parameterized simulations
\vfill
\hfill {\tiny CoDA 2014 {\em Statistical Analysis and Data Mining} tex template (do not remove)}

\newpage
\spacingset{1.45}

\section{Introduction}
\label{sec:intro}

Science and engineering simulations often contain several input parameters---e.g., physical constants, boundary conditions, or geometry descriptions. When presented with such parameterized simulations, the scientist naturally wonders how the simulated predictions depend on the input parameters. Which parameters, when perturbed, create the largest change in predictions? How precisely must the parameters be specified to ensure accurate predictions? And what is the effect of imprecisely prescribed input parameters on the predictions? If the goal is to maximize or minimize a predicted quantity, which combinations of parameters correspond to larger or smaller values of the prediction? In a similar vein, one may ask which parameter values correspond to predictions outside a safe region of operation---or which parameter values yield predictions that are consistent with a set of observations.

Many scientists rely on intuition about the physical system to answer these questions. But intuition becomes less trustworthy as the simulations become more complex, e.g., when they include several interacting physical components. Algorithms for optimization, uncertainty quantification, and model calibration become more attractive as model complexity increases. If an algorithm can easily interface with the simulation---e.g., by automatically evaluating predictions given values for the inputs---then applying the algorithm to the simulation becomes relatively easy.

The number of times the algorithm needs to evaluate a prediction increases---sometimes extremely rapidly---as the number of inputs increases. The situation is worse if each evaluation requires significant computational resources. For example, finding the global minimum of a complicated prediction depending on 100 inputs is not tractable if available resources permit only ten model evaluations. In practice, the scientist may choose only the most important parameters to vary so the study fits within the computational budget. Alternatively, one may seek a low-dimensional description of the prediction as a function of the input parameters. If such a description is sufficiently accurate, then studies like optimization or calibration can work in the space of fewer variables---potentially allowing the desired parameter studies within the given budget.

The \emph{active subspace} is a low-dimensional linear subspace defined by important directions in the model's input space; input perturbations along these directions change the model's prediction more, on average, than perturbations orthogonal to the important directions. Not all models have an active subspace. Some model's predictions respond significantly to input perturbations along all directions. However, if a model does admit an active subspace, then one can exploit it to perform parameter studies in the coordinates of the subspace---i.e., the \emph{active variables}---which are linear combinations of normalized versions of the model's input parameters. Therefore, it can be very advantageous to discover that a model admits an active subspace.

Active subspaces have been studied in a variety of contexts under different names. Cook's excellent text \emph{Regression Graphics}~\cite{cook2009regression} reviews and develops statistical methods for dimension reduction in the context of regression surfaces, and it contains references to the major works in the statistics literature. What we call the \emph{active subspace} is a type of \emph{dimension-reduction subspace} in Cook's parlance~\cite[Chapter 6]{cook2009regression}---though we are working with noiseless computer simulations in contrast to general regression surfaces. Russi's 2010 Ph.D.~dissertation uses the phrase \emph{active subspace} in a way comparable to our use~\cite{Russi2010}. He exploits the active subspace to construct quadratic surrogate models for uncertainty quantification in chemical kinetics models. Our prior work develops a theoretical framework for active subspaces including applications to kriging response surfaces~\cite{Constantine2013}. We have applied these methods to several models in aerospace engineering~\cite{constantine2011,wang2011,lukaczyk2014}. Abdel-Khalik has applied similar methods in nuclear engineering applications~\cite{AbdelKhalik2013}.

In this paper, we describe how to test a model for an active subspace, and we apply this test to the maximum power from a single-diode model of a photovoltaic solar cell with five input parameters. In Section \ref{sec:asm}, we generically describe the active subspace and how to search for it. We then describe the single-diode model, its input parameters, and its predicted \emph{performance parameters} (i.e., model outputs) in Section \ref{sec:pv}. We apply the tests for the active subspace to the model's maximum power in Section \ref{sec:pvasm} and show that a dominant one-dimensional active subspace is present in the five-dimensional space of input parameters. We conclude in Sections \ref{sec:discussion} and \ref{sec:conclusion} with a discussion of how one may exploit the low-dimension of the active subpace to further study the parameter dependence in the single-diode model's maximum power.

\section{Active subspaces}
\label{sec:asm}

We consider a generic multivariate function $f=f(\vx)$, where $\vx$ represents the inputs of the model, and $f$ represents a specific scalar performance parameter that the model predicts. Let $\sX=[-1,1]^m$ be the domain with $\vx\in\sX$, and let $\rho:\sX\rightarrow\mathbb{R}_+$ be a bounded and continuous weight function on $\sX$; we assume $\rho$ is normalized to integrate to 1. Assume $f$ is diffentiable and absolutely continuous, and denote the gradient $\nabla_\vx f(\vx) = [\partial f/\partial x_1,\dots,\partial f/\partial x_m ]^T$ oriented as a column vector.

Consider the following matrix $\mC$ defined as
\begin{equation}
\label{eq:C}
\mC \;=\; \int_\sX (\nabla_\vx f)(\nabla_\vx f)^T\,\rho\,d\vx.
\end{equation}
In the context of dimension reduction for regression functions, Samarov calls this matrix an \emph{average derivative functional}~\cite{samarov1993exploring}. Note that we are not studying regression functions, per se. Instead, $f$ represents a noiseless, parameterized computer simulation. The matrix $\mC$ is symmetric and positive semi-definite, so it admits a real eigenvalue decomposition
\begin{equation}
\mC \;=\; \mW\Lambda\mW^T,
\qquad
\Lambda \;=\; \mathrm{diag}(\lambda_1,\dots,\lambda_m),
\qquad
\lambda_1\geq\cdots \geq\lambda_m\geq 0.
\end{equation}
We can partition the eigenvectors,
\begin{equation}
\label{eq:partition}
\mW = \bmat{\mW_1 & \mW_2},\qquad \Lambda = \bmat{\Lambda_1 & \\ & \Lambda_2},
\end{equation}
where $\mW_1$ contains the first $n$ eigenvectors, and $\Lambda_1$ contains the $n$ largest eigenvalues. We use the two sets of eigenvectors to create new sets of variables $\vy=\mW_1^T\vx$ and $\vz=\mW_2^T\vx$. We call the subspace defined by $\mW_1$ the \emph{active subspace}, and we call the variables $\vy$ the \emph{active variables}; the term \emph{active subspace methods} for this type of analysis was first used in Russi's 2010 Ph.D.~dissertation~\cite{Russi2010}. The following two lemmas justify these labels.

\begin{lemma}
\label{lem:avgsqgrad}
The mean-squared directional derivative of $f$ with respect to the eigenvector $\vw_i$ is equal to the corresponding eigenvalue,
\begin{equation}
\int_\sX \big((\nabla_{\vx} f)^T\vw_i\big)^2\,\rho\,d\vx \;=\; \lambda_i.
\end{equation}
\end{lemma}

\begin{lemma}
\label{lem:grad}
The mean-squared gradients of $f$ with respect to the coordinates $\vy$ and $\vz$ satisfy
\begin{equation}
\begin{aligned}
\int_\sX (\nabla_{\vy} f)^T(\nabla_{\vy} f)\,\rho\,d\vx &= \lambda_1+\cdots+\lambda_n, \\
\int_\sX (\nabla_{\vz} f)^T(\nabla_{\vz} f)\,\rho\,d\vx &= \lambda_{n+1}+\cdots+\lambda_m.
\end{aligned}
\end{equation}
\end{lemma}
The proofs of these lemmas can be found in our prior work~\cite{Constantine2013}. In words, they mean that $f$ changes more, on average, when its inputs are perturbed along the directions $\mW_1$ than along the directions $\mW_2$; the eigenvalues quantify precisely how much more. If an eigenvalue is exactly zero, then $f(\vx)$ is constant along the direction defined by the corresponding eigenvector over all of $\sX$.

To gain some intuition, consider the extreme case where all eigenvalues are precisely zero except $\lambda_1$. Then $f(\vx) = g(\vw_1^T\vx)$, where $\vw_1$ is the first column of $\mW$, and $g$ is a function of one variable. In many applications, the smaller eigenvalues are not precisely zero, but they may be much (e.g., orders of magnitude) smaller so that $f(\vx)$ may be reasonably approximated by a function of $n<m$ linear combinations of $\vx$.

If a given model $f$ admits such structure, then certain operations---e.g., response surface modeling or optimization---become much less expensive. In particular, these operations can be performed in the $n$-dimensional space of the active variables $\vy$ instead of the full $m$-dimensional space. It is therefore extremely valuable to determine if $f$ admits an active subspace. This analysis assumes that the weight function $\rho$ on the domain $\sX$ is given. The computed quantities (like $\mW$ and $\Lambda$) will change if a different $\rho$ is given.

\subsection{Identifying an active subspace}
To identify the active subspace, we must approximate the matrix $\mC$ from \eqref{eq:C}. Since $\mC$ is the mean of the outer product of the gradient, we can approximate it with simple Monte Carlo. One could achieve a more accurate approximation of $\mC$ and its eigendecomposition with an integration rule that is more accurate than simple Monte Carlo. However, if $m$ is greater than two or three, then tensor product constructions of accurate univariate numerical integration rules (e.g., Gaussian quadrature) require too many evaluations of the gradient to be practical---especially if the gradient is expensive to compute. The simple Monte Carlo proceeds as follows. Draw $\vx_i$ independently according to $\rho$ with $i=1,\dots,M$. Then
\begin{equation}
\label{eq:Cmc}
\mC \;\approx\;\hmC\;=\;
\frac{1}{M}\sum_{i=1}^M \nabla_\vx f(\vx_i)\,\nabla_\vx f(\vx_i)^T \;=\;
\hmW\hLambda\hmW^T.
\end{equation}
We use the eigenvalues of the approximation $\hmC$ as evidence of an active subspace. In particular, a large gap in the eigenvalues indicates a separation between the corresponding active and inactive subspaces defined by $\hmW_1$ and $\hmW_2$, respectively. The qualification \emph{large} depends on the application. For example, a good low-dimensional approximation of $f$ over its entire domain may need a much larger separation than a good approximation of the bounds or the average of $f$.

Other eigenvalue-based dimension reduction schemes (e.g., principal component analysis) use heuristics based on the magnitude of the eigenvalues, such as choosing the dimension of the subspace such that the ratio $\sum_{i=1}^n \lambda_i/\sum_{i=1}^m \lambda_i$ is larger than some threshold. For active subspaces, such heuristics are not well-justified. We are not interested in capturing the gradient's variance; we want to identify directions along which the function changes the most. In our recent work, we show that the distance between the $n$-dimensional subspace defined by $\mW_1$ and its approximation defined by $\hmW_1$ is inversely proportional to $\lambda_{n+1}-\lambda_n$~\cite{constantine2014computing}. Therefore, the quality of the finite-sample subspace approximation depends more on the gap between the eigenvalues than their magnitude.

\subsection{Bootstrap to estimate variability}
\label{sec:boot}
If a gap is present in the approximate eigenvalues $\hLambda$, one may naturally ask if a comparable gap is present in the true eigenvalues $\Lambda$. To address this question, we use nonparametric bootstrap. The bootstrap is most appropriate when the simulation is expensive, and the number $M$ of gradient evaluations is constrained by a computational budget. For $j=1,\dots,M'$, let $\bpi_j = [\pi_1^j,\dots,\pi_M^j]$ be an $M$-vector of integers drawn uniformly at random between 1 and $M$. The $j$th \emph{bootstrap replicate} of $\hmC$ is computed as
\begin{equation}
\label{eq:Cboot}
\hmC^j\;=\;
\frac{1}{M}\sum_{i=1}^M \nabla_\vx f(\vx_{\pi_i^j})\,\nabla_\vx f(\vx_{\pi_i^j})^T \;=\;
\hmW^j\hLambda^j\left(\hmW^j\right)^T.
\end{equation}
The collection of eigenvectors $\hLambda_j$ yields a bootstrap distribution for the eigenvalues, which can be used to estimate boostrap intervals. Note that these intervals are not true confidence intervals, since the estimates are biased. However, this bias decreases as $M$ increases; see~\cite[Section 7.2]{efron1994introduction} and~\cite[Chapter 3]{jolliffe2002principal} for related discussions of the bias in bootstrap estimates for principal components. Nevertheless, we use the 99\% bootstrap intervals as evidence of gaps in the estimated eigenvalues.

We can also use the bootstrap replicates to estimate the error in the estimated subspace. Partition the eigenvectors from the bootstrap replicate as in \eqref{eq:partition}. Define
\begin{equation}
\label{eq:subspacerr}
e_j \;=\; \|\hmW_1\hmW_1^T - (\hmW^j_1)(\hmW^j_{1})^T\|_2
\;=\; \|\hmW_1^T\hmW^j_{2}\|,
\end{equation}
which is the distance between the subspace defined by $\hmW_1$ and bootstrap replicate $\hmW_1^j$~\cite{golub1996matrix}. The mean and 99\% bootstrap intervals of the set $\{e_j\}$ quantify the varibility in the estimated active subspace.

Our current research efforts are focused on rigorously justified, computable metrics for determining the relationship between the approximate eigenpairs $\hLambda$, $\hmW$ and the true eigenpairs $\Lambda$, $\mW$. Our recent paper applies non-asymptotic random matrix theory to characterize the error in these approximations~\cite{constantine2014computing}.

\subsection{Visualization with the active subspace}
Scatter plots are a common way to visualize data sets in search of a trend. Unfortunately, visualization tools can only display scalar responses $f$ as a function of at most two variables. When the response depends on more than two variables, one can plot responses versus each variable or each pair of variables.

Using the vectors defining the active subspace, we can create scatter plots based on the active variables, which often show discernible trends that can be exploited in further studies (e.g., building response surfaces or optimization). This idea is described in detail in the texts on \emph{regression graphics}~\cite{cook2009regression}, where such plots are called \emph{sufficient summary plots}. For example, suppose we have noticed a large gap between the first and second eigenvalue. We can confirm the presence of the active subspace by first sampling $\vx_i$ from $\rho$ as in \eqref{eq:Cmc}, computing $f_i = f(\vx_i)$ (often computed along with the gradients in \eqref{eq:Cmc}), and plotting the pairs $(\vw_1^T\vx_i,f_i)$---where $\vw_1$ is the first column of $\mW$. If a tight, univariate trend is visibly present, then this verifies the active subspace. Figure \ref{fig:proj1} in Section \ref{sec:pvasm} shows these scatter plots for the single-diode model.

These sufficient summary plots can also be used to examine the variability in the vectors defining the active subspace. In particular, we can use the bootstrap replicates $\hmW^j$ from \eqref{eq:Cboot}. The first eigenvector from each can be used to create pairs for the scatter plot. Plotting all pairs together creates clusters in the scatter plot that can indicate the spread around a trend. We create these visualizations for the predicted performance parameter from the single-diode model in Section \ref{sec:pvasm}.

\subsection{Gradient approximation}
Often the gradient $\nabla_\vx f$ is not available or is too complicated to compute. To approximate partial derivatives, one builds a model of $f(\vx)$ that is easily differentiable, e.g., a polynomial model. A first-order finite difference approximation of the partial derivative with respect to $x_i$ at $\vx$ computes the slope of the plane that interpolates $f$ at $\vx$ and $\vx+\varepsilon \ve_i$, where $\ve_i$ is a vector of zeros with a one in the $i$th component. When the variables are properly normalized, a reasonable choice for $\varepsilon$ to approximate the derivative and avoid round-off issues is the square root of machine precision~\cite[Section 8.1]{Nocedal2006}.

If the simulation output is well-behaved as a function of the variables $\vx$, then finite difference approximations of gradients can be used as a substitute for the true gradient. Each finite difference approximation requires $m+1$ evaluations of $f$, so the number of evaluations to compute $\hmC$ is $M(m+1)$. Therefore, finite differences are most appropriate when the following are satisfied:
\begin{itemize}
\item the simulation output behaves well as a function of the parameters, i.e., there is no noise due to limited iterations of a nonlinear solver or large changes of the function on scales smaller than the finite difference step size,
\item the simulation completes in a short enough time on available computing resources to permit $M(m+1)$ evaluations.
\end{itemize}
Both of these conditions are satisfied for the performance parameter of interest from the single-diode model as a function of its input parameters. Therefore, we use finite differences in place of the true gradient to search for an active subspace.

If these conditions are not satisfied, then one must pursue other options for estimating the gradients. One idea is to use local linear models within a budget-constrained set of runs. In short, assume one has $B$ pairs $(\vx_j,f(\vx_j))$. Given a point $\vx$ in the parameter space, select $M$ points in $\{\vx_j\}$ near $\vx$, fit a linear model to the selected pairs, and estimate the gradient as the gradient of the linear model. $M$ must be large enough and the points must be chosen to enable fitting the linear model. We have used this approach in practice, but we are still analyzing its properties. The advantage over finite differences is that it is robust to noise in $\{f_j\}$ when the number of points used to fit the linear model is greater than the number of variables ($M>m+1$). However, if the points $\{\vx_j\}$ are not sufficiently dense in the parameter space, then the gradient estimates can be very poor, which lead to poor estimates of the eigenvalues and eigenvectors. We continue to analyze this approach.

\section{A single-diode solar cell model}
\label{sec:pv}

In this section we introduce the single-diode solar cell model, the input parameters, and the predicted quantity of interest. Over the past decade, a growing demand for clean energy has caused rapid growth in the photovoltaic (PV) industry. The overall health of the industry depends on proper characterization of the risks associated with a PV system. One can characterize such risks using a mathematical model to describe the performance of the PV device, such as a single-diode lumped-parameter equivalent-circuit model defined in Section \ref{sec:params}. For many series-wired PV devices, the single-diode model accurately describes the device's current-voltage ($I$-$V$) characteristics at given irradiance and temperature conditions. The single-diode model contains several input parameters that must be estimated---e.g., with methods in \cite{Hansen2013b} or \cite{Farivar2010}---to properly pose the model. With the inputs fixed, practitioners can estimate the {\it key performance parameters}, such as the maximum power output $\Pmax$ and energy conversion efficiency $\eta$.

Estimating the model's input parameters presents challenges. For example, many single-diode model parameter estimation methods have been proposed in the PV literature, and it is not always clear if any particular estimation yields a unique result \cite{Hansen2013}. Further, since each key performance parameter is a function of the model parameters, errors in estimating the latter will cause issues in estimating the former.

Given these issues, it is desirable to know how the variability or uncertainty in each of the input parameters affects variability in the outputs. In particular, it would be useful to know whether the model admits an active subspace. If an active subspace is present, then estimation of the performance parameters is greatly simplified. For example, rather than conducting current-voltage measurements in a way that yields accurate estimates of all model parameters, practitioners could instead focus on measurements that accurately estimate parameters in the active subspace.

\subsection{Calculating the key performance parameters}
\label{sec:params}

The single-diode model describes a relationship between the current ($I$) and voltage ($V$) in the single diode. There are multiple variants of the single-diode model, each of which has accompanying auxiliary equations that describe how components depend on irradiance and temperature. In this work, we use the auxiliary equations derived in \cite{Zaharatos2014}; see \cite{DeSoto2006} for an example of others. These relationships are defined in the following equations.
\begin{equation}
\label{eq:model}
I  \;=\; \IL - \IS\,\left(\exp\left(\frac{V + I\,\RS}{\NS\,n \,\Vth}\right) - 1 \right) - \frac{V + I\,\RS}{\RP}.
\end{equation}
$\NS$ is the number of cells connected in series, which we set to 1. The thermal voltage $\Vth$ is a given, fixed constant for fixed temperature $T=25^\circ\mathrm{C}$. An auxiliary equation defines the photocurrent $I_\text{L}$ as
\begin{equation}
\label{eq:photocurrent}
I_\text{L} \;=\; \ISC + \IS \left(\exp\left(\frac{\ISC \RS}{\NS\, n\, \Vth}\right) - 1 \right) + \frac{\ISC\,\RS}{\RP}.
\end{equation}
The remaining terms are input parameters to the single-diode model; their names, units, and ranges are given in Table \ref{tab:domains}. The ranges correspond to typical 2cm$^2$ crystalline silicon PV cells.

\begin{table}[ht]\small
\caption{The ranges of the input parameters for the single-diode model in \eqref{eq:model} and \eqref{eq:photocurrent}. These ranges correspond to typical 2cm$^2$ crystalline silicon PV cells.}
\begin{center}
\begin{tabular}{lllll}
Parameter & Name & Lower bound & Upper bound & Units\\
\hline
$\ISC$ & short-circuit current & 0.05989 & 0.23958 & amps\\
$\IS$ & diode reverse saturation current & 2.2e-11 & 2.2e-7 & amps \\
$n$ & ideality factor & 1 & 2 & unit-less\\
$\RS$ & series resistance & 0.16625 &  0.66500 & ohms \\
$\RP$ & parallel (shunt) resistance & 93.75 & 375.00 & ohms
\end{tabular}
\end{center}
\label{tab:domains}
\end{table}

Given a device and $I$-$V$ data, one can estimate precise values for the input parameters $(\ISC, \IS, n, \RS, \RP)$; see \cite{Hansen2013b,Farivar2010,Zaharatos2014}. With these values fixed, one can compute the maximum power of the device as
\begin{equation}
\label{eq:power}
\Pmax \;=\; \maxi{I,V}\;I\, V,
\end{equation}
where current $I$ and voltage $V$ are constrained by \eqref{eq:model} and \eqref{eq:photocurrent}. Changing the inputs $(\ISC, \IS, n, \RS, \RP)$ changes the nonlinear constraints in the optimization \eqref{eq:power}, thus changing $\Pmax$. Therefore, we can write
\begin{equation}
\Pmax \;=\; \Pmax(\ISC, \IS, n, \RS, \RP).
\end{equation}
We use MATLAB to solve the optimization \eqref{eq:power} subject to the constraints \eqref{eq:model} and \eqref{eq:photocurrent}. Figure \ref{fig:ivcurves} shows ten $I$-$V$ curves with the single-diode model parameters chosen uniformly at random from the ranges in Table \ref{tab:domains}. The black dots show the values of $I_{\text{max}}$ and $V_{\text{max}}$ that produce $\Pmax$. A wrapper function computes $\Pmax$ as a function of the input parameters. Requests for the MATLAB code may be sent to the third author Campanelli.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.45\textwidth]{ivcurves-crop.pdf}
\end{center}
\caption{Ten current-voltage curves with parameters chosen uniformly at random from the ranges in Table \ref{tab:domains}. The black circles plot the location of $\Pmax=I_{\text{max}}V_{\text{max}}$.}
\label{fig:ivcurves}
\end{figure}

\section{An active subspace in the single-diode model}
\label{sec:pvasm}

We apply the techniques from Section \ref{sec:asm} to search for an active subspace of the $\Pmax$ performance parameter as a function of the parameters defined in Table \ref{tab:domains}. Note that a different performance parameter (e.g., efficiency) would require an independent and identical analysis. The relationship between the respective active subspaces depends on the relationship between the performance parameters and their inputs. In some cases, the active subspaces might be similar; in other cases, they may differ substantially.

The MATLAB scripts for generating the figures in this section can be found at \url{https://bitbucket.org/paulcon/active-subspaces-in-a-single-diode-model}. To generate the figures and results, we evaluated the model $(m+1)M=6000$ times, which took approximately three minutes on a MacBook Air. A more expensive model would permit fewer evaluations, and we would choose $M$ differently---or perhaps choose a different method.

\subsection{Normalizing the input parameters}

We first normalize the domain of $\Pmax$ to the hypercube; denote the normalized variables by $\vx=[x_1,x_2,x_3,x_4,x_5]^T$. The variable $\IS$ varies over several orders of magnitude, and preliminary tests show that $\Pmax$ changes rapidly near smaller values of $\IS$. To address this, we work with $\log(\IS)$. The range of $\log(\IS)$ is bounded below by -24.54 and above by -15.32.

Let $\sL(\cdot)$ and $\sU(\cdot)$ return the upper and lower bounds, respectively, of the argument. We define the normalized input parameters as
\begin{equation}
\begin{aligned}
x_1 &= 2\,\left( \frac{\ISC-\sL(\ISC)}{\sU(\ISC)-\sL(\ISC)} \right) - 1\\
x_2 &= 2\,\left( \frac{\log(\IS)-\sL(\log(\IS))}{\sU(\log(\IS))-\sL(\log(\IS))} \right) - 1\\
x_3 &= 2\,\left( \frac{n-\sL(n)}{\sU(n)-\sL(n)} \right) - 1\\
x_4 &= 2\,\left( \frac{\RS-\sL(\RS)}{\sU(\RS)-\sL(\RS)} \right) - 1\\
x_5 &= 2\,\left( \frac{\RP-\sL(\RP)}{\sU(\RP)-\sL(\RP)} \right) - 1
\end{aligned}
\end{equation}
Then $\vx\in[-1,1]^5$ and has no units. Note that this map is invertible. In other words, given a value $\vx\in[-1,1]^5$, one can shift and scale the components of $\vx$---and transform the second component with the exponential---to produce a valid input for the single-diode model. We take the weight function $\rho$ in \eqref{eq:C} to be a constant $2^{-5}$ inside $[-1,1]$ and zero elsewhere. One can interpret this weight function as a uniform probability density on the space of input parameters.

Proper scaling is important. The results of the active subspaces analysis are not scale invariant, so the scientist must carefully choose input parameter bounds that are appropriate for the application. Changing the bounds can change the results---sometimes dramatically. For example, the quantity of interest may be irregular and badly behaved over a wide range of inputs. But reducing the range of interest might focus on a small region in the parameter space where $f$ is smooth---even nearly linear. A linear function has a one-dimensional active subspace. Thus, loosely speaking, if the range of parameters is sufficiently small, then the active subspace will be nearly one-dimensional for continuous and differentiable quantities of interest.

\subsection{Estimating the eigenvalues and eigenvectors}
We choose $M=1000$ points $\vx_i$ uniformly at random from $[-1,1]^5$, and for each $\vx_i$ we compute both $\Pmax$ and a first-order finite difference approximation of the gradient. We use a finite difference step size of $10^{-6}$ in the normalized domain, which was small enough to ensure that gradient approximation errors were substantially smaller than the variance due to random sampling.

To study variability in the computed components of $\hmC$ from \eqref{eq:Cmc}, we use the bootstrap described in Section \ref{sec:boot} with $M'=1000$ bootstrap replicates. The 99\% boostrap intervals for the eigenvalue estimates are shown in Figure \ref{fig:evals}. The tight bootstrap ranges suggest there is little eigenvalue variability in the bootstrapped data set. There is a gap of nearly an order of magnitude between both the first and second and the second and third eigenvalues. This suggests a very dominant one-dimensional active subspace. Figure \ref{fig:se} shows the estimated error in the approximated subspace using \eqref{eq:subspacerr}. Note that this error is actually the cosine of the principal angle between the subspaces, so it is bounded above by 1.

\begin{figure}[h]
\begin{center}
\subfloat[Eigenvalues]{\label{fig:evals}
\includegraphics[width=0.45\textwidth]{evals-crop.pdf}
}
\subfloat[Subspace error]{\label{fig:se}
\includegraphics[width=0.45\textwidth]{se-crop.pdf}
}
\end{center}
\caption{The left figure shows the eigenvalues of the matrix $\hmC$ from \eqref{eq:Cmc} with $M=1000$ samples and their bootstrap ranges with $M'=1000$ boostrap replicates. The right figure shows the mean and ranges of the bootstrap estimates of the subspace error (see \eqref{eq:subspacerr}) with $M'=1000$ bootstrap replicates.}
 \label{fig:evals-and-se}
\end{figure}

Figure \ref{fig:evecs} displays the components of the first and second eigenvectors of $\hmC$ computed with $M=1000$ samples---along with the $M'=1000$ bootstrap replicates. They are normalized so that the first component is positive. The small ranges suggest that the one- and two-dimensional active subspaces are stable within the gradient samples.

\begin{figure}[h]
\begin{center}
\subfloat[First eigenvector]{
\includegraphics[width=0.45\textwidth]{w1-crop.pdf}
}
\subfloat[Second eigenvector]{
\includegraphics[width=0.45\textwidth]{w2-crop.pdf}
}
\end{center}
\caption{These figures show the components of the first (left) and second (right) eigenvectors of $\hmC$ with $M=1000$ gradient samples. The grey shaded region superimposes the $M'=1000$ bootstrap replicates of the eigenvector components; the tight ranges suggests stability in the gradient samples. These values can be used as measures of sensitivity of $\Pmax$ with respect to the input variables in the model, whose names are given along the horizontal axis of the figures. See Table \ref{tab:domains} for a description of the variables.}
 \label{fig:evecs}
\end{figure}

\subsection{Identifying the active subspace}
The order-of-magnitude gap between the first and second eigenvalues suggests a dominant one-dimensional active subspace. Recall that gaps between eigenvalues are more important than the magnitudes of the eigenvalues. Figure \ref{fig:oned} plots 100 values of $\Pmax$ against a linear combination of the corresponding normalized inputs; the weights of the linear combination are the components of the first eigenvector $\hvw_1$ from \eqref{eq:Cmc}. Figure \ref{fig:oned} is a one-dimensional sufficient summary plot, as described in~\cite{cook2009regression}. The plot shows a potentially exploitable low-dimensional relationship between the input parameters and the quantity of interest, $\Pmax$.

\begin{figure}[h]
\begin{center}
\subfloat[Projected with first eigenvector]{
\includegraphics[width=0.45\textwidth]{oned-crop.pdf}
\label{fig:oned}
}
\subfloat[Projected with first two eigenvectors]{
\includegraphics[width=0.45\textwidth]{twod-crop.pdf}
\label{fig:twod}
}
\end{center}
\caption{One-dimensional (left) and two-dimensional sufficient summary plots of 100 realizations of $\Pmax$ using $\hvw_1$ (left) and both $\hvw_1$, $\hvw_2$ (right). The clusters of points on the right correspond to 20 bootstrap replicates of $\hvw_1$ and $\hvw_2$. The strong trend in \ref{fig:oned} is verified by looking at the variation along \emph{Active Variable 1} (the horizontal axis) in \ref{fig:twod}, and the combination of these provide evidence of a one-dimensional active subspace.}
\label{fig:proj1}
\end{figure}

Figure \ref{fig:twod} shows the two-dimensional sufficient summary plot. The grayscale corresponds to the value of $\Pmax$. The black dot in the center of each cluster has a horizontal and vertical location determined by linear combinations of $\Pmax$'s corresponding input parameters; the weights of the linear combination are the eigenvectors $\hvw_1$ and $\hvw_2$, respectively. Each cluster corresponds to the value of $\Pmax$ with horizontal and vertical position determined by 20 bootstrap replicates of $\hvw_1$ and $\hvw_2$. The relatively small spread in the clusters indicates that the estimated subspace is stable.

In Figure \ref{fig:twod}, note that $\Pmax$ increases primarily as a function of the first active variable. This provides more evidence of the dominance of the one-dimensional active subspace.

\subsection{Eigenvectors for sensitivity analysis}

The eigenvector components connect the active subspace to the (normalized) variables in the model; Figure \ref{fig:evecs} includes the corresponding labels in the horizontal axis. The magnitudes of the eigenvector components can be used as measures of relative sensitivity for each of the parameters in the model. A large absolute value of the eigenvector component implies that this variable is important in defining the direction along which input perturbations produce the most change, on average, in $\Pmax$. The components corresponding the parallel resistance $\RP$ and series resistance $\RS$ are close to zero, which implies that (normalized) changes in $\RP$ and $\RS$ do not change $\Pmax$ as much as changes in the other parameters.

We compare the components of $\hvw_1$, which defines the dominant one-dimensional active subspace, to the Sobol' first-order sensitivity indices and total sensitivity indices. We implemented the method to estimate these indices from Sudret~\cite{Sudret2008} using a tensor product Gaussian quadrature rule with 32768 points in 5 dimensions. Table \ref{tab:sens} displays the eigenvector components and the Sobol' indices.

\begin{table}[ht]
\caption{The components of the eigenvector $\hvw_1$, the Sobol' first-order indices, and the Sobol' total sensitivity indices for the single-diode model parameters. The metrics provide same ranking of importance for the parameters.}
\begin{center}
\begin{tabular}{cccc}
Parameter & $\hvw_1$ & Sobol' first-order index & Sobol' total sens.~index\\
\hline
$\ISC$ & 0.77 & 0.56 & 0.61\\
$\IS$ & -0.42 & 0.17 & 0.19\\
$n$ & 0.47 & 0.21 & 0.25\\
$\RS$ & -0.09 & 0.01 & 0.01\\
$\RP$ & 0.02 & 0.00 & 0.00
\end{tabular}
\end{center}
\label{tab:sens}
\end{table}

There are several things to note in Table \ref{tab:sens}. First, the importance ranking from the Sobol' indices and the ranking from the eigenvector component magnitudes are the same. However, the signs of eigenvector components indicate whether $\Pmax$ will increase or decrease, on average, with changes in the corresponding parameter. (Eigenvectors are unique up to a sign, so these signs should be considered relative to one another.) There is no such interpretation in the Sobol' indices. In other words, the eigenvector components provide more information about the relationship between the corresponding parameters and $\Pmax$.

The numeric values are difficult to compare, since they are normalized differently. The Sobol' indices are divided by the estimated total variance in $\Pmax$, whereas the eigenvector components are normalized to have Euclidean norm equal to 1. In this case, they are the same order of magnitude.

We emphasize the different interpretations of these numbers. The eigenvector $\hvw_1$ is a single direction in the input space. Input perturbations along this direction change $\Pmax$ more, on average, than perturbations orthogonal to this direction. The Sobol' indices indicate the proportion of $\Pmax$'s total variance attributable to the factors of a variance-based decomposition of $\Pmax$. In terms of dimension reduction, one might (i) use the Sobol' indices to conclude that only three of the five parameters were important and (ii) approximate $\Pmax$ as a function of those three parameters. In contrast, the active subspace is used to approximate $\Pmax$ as with a univariate function of a linear combination of all five parameters. In other words, the variance-based approach reduces the dimension from 5 to 3, and the active subspace approach reduces the dimension from 5 to 1.

There are other challenges facing the active subspace style of dimension  that do not affect the variance-based approach. The variance-based approach uses a subset of the model's parameters---as opposed to a subspace. Fixing a subset of parameters at nominal values and allowing the others to vary is relatively straighforward. Exploiting one important linear combination of the normalized parameters is not as straightforward. In the next section we discuss some possibilities for exploiting the active subspace.

\subsection{Interpreting the sensitivity analyses}
We can interpret the consistent importance rankings in the sensitivity metrics in terms of the physics of the single-diode model. First, $\Pmax= I_{\text{max}}V_{\text{max}}$, where $(I_{\text{max}},V_{\text{max}})$ is the point on the $I$-$V$ curve that maximizes \eqref{eq:power}. $I_{\text{max}}$ typically scales linearly with $\ISC$. The voltage at $I=0$---denoted $V_{\text{oc}}$ for \emph{voltage at open circuit}---typically scales logarithmically with $\ISC$. $V_{\text{max}}$ behaves like $V_{\text{oc}}$. Thus, increases in $\ISC$ affect $I_{\text{max}}$ linearly and $V_{\text{max}}$ logarithmically, which changes $\Pmax$ as the product of a linear term and a logarithmic term. Since $\ISC$ affects both $I_{\text{max}}$ and $V_{\text{max}}$, it is natural for it to be the driving parameter in $\Pmax$.

Recall that these sensitivity analyses are not scale invariant; the input parameter ranges and the weight function affect the results. The ranges for $n$ and $\IS$ are at about their maximum extent for a crystalline silicon PV cell, while the ranges for $\RS$ and $\RP$ are somewhat restricted. All of these parameters affect the shape and/or onset of the \emph{knee} of the $I$-$V$ curve, which determines the location of $(V_{\text{max}}, I_{\text{max}})$ and, thus, $\Pmax$. Here $\RS$ and $\RP$ are relatively unimportant because their ranges are relatively smaller than the ranges for $n$ and $\IS$.

\section{Using the active variables}
\label{sec:discussion}

The estimated eigenvalues in Figure \ref{fig:evals} and the plots in Figures \ref{fig:oned} and \ref{fig:twod} provide strong evidence for the presence of a dominant one-dimensional active subspace in the $\Pmax$ performance parameter computed from the single-diode model as a function of the model's five input parameters. The natural question is how one can exploit the dimension reduction afforded by the active subspace. There are five types of studies that benefit greatly from fewer input parameters. To keep the scope of this paper limited, we do not address the details of any of these studies, and we prefer to reserve them for future exploration.

\noindent \textbf{Visualization.} The active subspace enables one to view the model output's dependence on its inputs with standard computer graphics tools when the active subspace is not more than two-dimensional. Sufficient summary plots, such as those in Figures \ref{fig:oned} and \ref{fig:twod} can provide insights to modelers seeking to improve their models.

\noindent \textbf{Optimization.} Suppose one wanted to maximize $\Pmax$ over the input variables. The plot in Figure \ref{fig:proj1} shows a monotonic trend in the univariate function of the first active variable. For such functions, maximization is trivial; simply make the first active variable $\hvw_1^T\vx$ as large as possible subject to the constraints $\vx\in [-1,1]^m$. In general, a five-dimensional global optimization where the input/output relationship is not well-understood is a difficult problem. For the $\Pmax$ output from the single-diode model, the active subspace provides a way to discover the location of the global optimum with ease.

\noindent \textbf{Response surfaces.} If the single-diode model were expensive to evaluate, then one may wish to construct a response surface that approximates the map from inputs to outputs. A one-dimensional active subspace allows one to build a response surface on only the active variable instead of the five model input parameters. Constructing a response surface in one variable is certainly preferred to constructing one in five variables. Real-time control systems such as maximum power point trackers might benefit from the characterization of a simple, low-dimensional response surface. In addition, the visualization tools give one confidence that the response is sufficiently smooth with respect to the active variable to permit an accurate response surface.

\noindent \textbf{Averages.} One may wish to compute an average of $\Pmax$ over all five input variables. Since the average of $\Pmax$ over its inputs is the same as the average of the conditional expectation given the active variable, we need only to approximate a marginal density of the active variable to compute the average. We are working on the specifics of this idea, but the goal would be to transform a five-dimensional integral into a one-dimensional integral requiring many fewer evaluations of $\Pmax$.

\noindent \textbf{Design.} Suppose a modeler wanted to design a solar cell with $\Pmax$ in a specified range and sought an appropriate range of values for the model inputs. A one-dimensional active subspace makes this query much easier---especially if $\Pmax$ is monotonic with respect to the active variable as in Figure \ref{fig:twod}.

\section{Conclusions}
\label{sec:conclusion}

We have discussed methods for discovering the directions in a model's input space that change the model's prediction the most, on average. The active subspace is the span of these directions. For global sensitivity analysis, the active subspace offers an alternative to variance decomposition techniques such as Sobol' indices. We have applied this procedure to a single-diode solar cell model with key performance parameter $\Pmax$, and we discovered a dominant one-dimensional active subspace. We offered several possible ways to exploit the knowledge of this low-dimensional parametric dependence to gain greater insight into the model.

\pagebreak

\bibliographystyle{natbib}
\bibliography{pv-asm-coda}

\end{document}


text:
\section{WFFs of PropLog of Length 1}
Tags: Language of Propositional Logic

\begin{theorem}
The only [[Definition:WFF of Propositional Logic|WFFs of propositional logic]] of [[Definition:Length of String|length]] $1$ are:
* The [[Definition:Letter|letters]] of the [[Definition:Formal Grammar of Propositional Logic|formal grammar of propositional logic]] $\mathcal L_0$
* The [[Definition:Top (Logic)|tautology symbol]] $\top$
* The [[Definition:Bottom (Logic)|contradiction symbol]] $\bot$.
\end{theorem}

\begin{proof}
We refer to the [[Definition:Bottom-Up Specification of Propositional Logic|rules of formation]].
From $\mathbf W: TF$, $\top$ and $\bot$ (both of [[Definition:Length of String|length]] 1) are WFFs.
From $\mathbf W: \mathcal P_0$, all elements of $\mathcal P_0$ (all of [[Definition:Length of String|length]] 1) are WFFs.
Every other rule of formation of the [[Definition:Formal Grammar of Propositional Logic|formal grammar of propositional logic]] consists of an existing WFF in addition to at least one other primitive symbol.
Hence the result.
{{qed}}
\end{proof}

config: wiki


text:
\begin{document}

\title{BSDEs driven by $G$-Brownian motion with  uniformly continuous generators}
\author{
Falei Wang \thanks{Zhongtai Securities Institute for Finance Studies and Institute for Advanced Research, Shandong University,
flwang@sdu.edu.cn. Research partially supported by  the National Natural Science Foundation of China  (No. 11601282) and the Natural Science Foundation of  Shandong Province (No. ZR2016AQ10).}
\and  Guoqiang Zheng\thanks{
School of Mathematics, Southeast  University, gqzheng@seu.edu.cn. Research partially  supported by the Fundamental Research Funds for the Central Universities (No. 3207018202) and
the National Natural Science Foundation of China (No. 11671231 and No. L162400033). } }
\maketitle
\date{}

\begin{abstract}
The present paper is devoted to investigating the existence and uniqueness of solutions to a class of non-Lipschitz  scalar valued  backward stochastic differential equations driven by $G$-Brownian motion ($G$-BSDEs). In fact, when the generators are Lipschitz continuous in $y$ and uniformly continuous in $z$, we construct the unique solution  to such equations by monotone convergence argument. The comparison theorem and related Feynman-Kac formula are stated as well.
\end{abstract}

\textbf{Keywords}: $G$-Brownian motion,  BSDE, uniformly continuous generators.

\textbf{Mathematics Subject Classification(2000)}: 60H10, 60H30.
\section{Introduction}
Given a Wiener space $(\Omega,\mathcal{F},\{\mathcal{F}_t\}_{t\in[0,\infty)},P_0)$, under which the canonical process $W_t$ constitutes a Brownian motion.
A typical nonlinear Lipschitz backward stochastic differential equations (BSDEs), which is formulated in Pardoux and Peng \cite{PP}, takes the form,
$$  Y_t=\xi +\int_t^T f(s,Y_s,Z_s)ds-\int_t^T  Z_s dW_s.$$
The authors found  a unique pair of adapted processes $(Y_s,Z_s)_{0\leq s\leq T}$ that satisfy the above equation for given squarely integrable  terminal value $\xi$ and Lipschitz generator $f$.

From then on, extensive  efforts have been made towards relaxing the Lipschitz  assumptions on the generator $f$. To mention just a few, for the scalar case, i.e., when $Y$ is 1-dimensional, Lepeltier and Martin \cite{LM} confirmed the existence of  solutions to BSDEs with continuous generator that is of linear growth. Kobylanski \cite{K} developed  the existence  for BSDEs with continuous generator that has a quadratic growth in $z$ when the terminal value $\xi$ is bounded. Also for the quadratic cases, Briand and Hu \cite{BH2006,BH2008} successively obtained the existence of solution for unbounded $\xi$.
More imaginative works on generalizing  the classical BSDEs theory from    different points of view   are emerging, and it would be too ambitious for us to give  an overview of all variants.
In this paper, we  focus on the study  of BSDEs under $G$-expectation framework.

The $G$-expectation theory was put forth by Peng \cite{P07,P08,P10}, which provides a unified  tool for stochastic analysis problems that involve non-dominated family of probability measures. In particular, the $G$-Brownian motion process is constructed with uncertain  quadratic variation process, a  feature that is helpful in capturing  the volatility fluctuations of financial market.
However there are many challenges in the research of $G$-expectation due to the uncertainty, for instance, the general dominated convergence theorem  is not available, see Example 11 in \cite{HJ1}.
Furthermore,  there exist non-increasing and continuous $G$-martingales called non-symmetrical $G$-martingales, which
 makes the $G$-martingale representation theorem  more difficult, see \cite{P10, STZ, Song11}.

 Recently, Hu, Ji, Peng and Song  \cite{HJPS2014} considered  the well-posedness problem of BSDEs driven by $G$-Brownian motion $B$:
 \begin{equation}\label{gbsd}
 Y_t=\xi +\int_t^T f(s,Y_s,Z_s)ds+\int_{t}^{T}g^{ij}(s,Y_s,Z_s)d\langle B^i, B^j\rangle_{s}-\int_t^T  Z_s dB_s- (K_T-K_t),
 \end{equation}
 where $\langle B^i, B^j\rangle$ denotes the quadratic (co)variation process and the generators $f,g$ are  Lipschitz continuous in $(y,z)$.
 The solution of  $G$-BSDE \eqref{gbsd} consists of a triple of processes $(Y,Z,K)$, where $K$ is a   non-symmetrical $G$-martingale.
 Note that  the classical Banach contraction mapping principle cannot be applied directly to this  equation due to the existence of $K$. The authors use PDE techniques and an   approximation of Galerkin type
 to obtain the existence and uniqueness result of $G$-BSDE \eqref{gbsd}. In an  accompany paper \cite{HJPS2014C}, they established the comparison theorem, Girsanov theorem and the relevant  nonlinear Feynman-Kac formula.

Note that there are at least two characteristics that make  the study  on $G$-BSDEs meaningful, for one thing, we could establish a connection with fully nonlinear parabolic partial differential equations (PDEs) using $G$-BSDEs, for the other, since there exists a family of non-dominated, mutually singular martingale measures underlying the $G$-Brwonian motion,  one can solve simultaneously a family of  classical BSDEs  driven by mutually singular continuous martingales through dealing with only one aggregated $G$-BSDE.  Moreover,  Song \cite{Song14} obtained gradient estimates for certain   nonlinear  partial differential equations (PDEs) by combining $G$-expectation theory with  coupling methods.
A close  approach to $G$-BSDEs is the so-called second order BSDEs framework proposed   independently by Soner, Touzi and Zhang \cite{STZ1}.

 Still there are further research papers on getting rid of the Lipschitz  assumptions, and extensions of $G$-BSDEs from different aspects. For instance, Hu, Lin and Soumana Hima \cite{HLH}
studied the $G$-BSDEs under quadratic assumptions on coefficients and Li, Peng and Soumana Hima \cite{LPH} considered $G$-BSDE with reflection,  for which situation the solution is forced
to lie above a prescribed continuous process. This paper is devoted to the research of the existence of solution to equation \eqref{gbsd} when  $f,g$ are Lipschitz continuous in $y$ and uniformly continuous in $z$, yet with a linear growth in both arguments.

In classical situation, the  term $B$ in \eqref{gbsd} boils down to the standard Wiener process, $d\langle B^i,B^j\rangle_t=(dt)\mathbf{1}_{i=j}$, and the term $K$ vanishes.
It was  due to Lepeltier and Martin \cite{LM} that confirmed those  equations  allow for solutions by monotone convergence argument, indeed their results hold for BSDEs that have continuous coefficients with linear growth,  and Jia \cite{J,J2010} supplemented with  proofs on uniqueness of solution.
However Lepeltier and Martin's arguments  cannot be applied directly to investigate the existence of solution to equation \eqref{gbsd},
 because the monotone convergence theorem  of $G$-expectation is hardly at hand and the convergence of approximating sequences of $G$-BSDEs is not obvious.

Our observation is,  this obstacle can be overcome with the help  of a uniform estimate for approximating sequences of $G$-BSDEs, see Lemma \ref{lemma2}.
Indeed, the  uniformly continuous generators can be approximated uniformly by a sequence of Lipschitz generators (see \cite{J2010}),
 from which we could prove the convergence of  approximating sequences of $G$-BSDEs based on the linearization method of \cite{HJPS2014C} and \cite{HW}.
Then we obtain  the existence and uniqueness of the solution to $G$-BSDE \eqref{gbsd} by $G$-stochastic analysis technique.
Since our work relies heavily on the comparison theorem of  Lipschitz $G$-BSDEs, we can only deal with the one dimensional case.
And the comparison theorem still holds for this type of $G$-BSDEs.
 Finally the connection of this equations with the second order fully nonlinear PDEs  are discussed,
 thanks to the stability of viscosity solution, we show the solution to Markovian $G$-BSDE \eqref{gbsd} defines the unique solution to the related  PDE, in the spirit of Feynman-Kac formula.

 The paper is organized as follows.  In the Section 2, we provide with preliminary notions on $G$-expectation and Lipschitz $G$-BSDEs.
 In Section 3, we state and prove our main theorem and the comparison theorem of our version. As an application, a slightly more general form of nonlinear Feynman-Kac formula is obtained in Section 4.

\section{Preliminaries}
To begin with, we shall recall some ingredients of $G$-expectation theory mainly from the seminal work of Peng \cite{P10}, and then of $G$-BSDEs results from \cite{HJPS2014,HJPS2014C}.

\subsection{$G$-expectation }
Consider the canonical path space $\Omega=C_{0}([0,\infty),\mathbb{R}^d)$, all continuous paths $\omega$ vanishing at zero, i.e., $\omega_0=0$. $(\Omega,\rho)$ is readily seen to be a complete separable metric space, where $\rho$ is given by,
\[
\rho(\omega^{1},\omega^{2}):=\sum_{i=1}^{\infty}2^{-i}[(\max_{t\in
\lbrack 0,i]}|\omega_{t}^{1}-\omega_{t}^{2}|)\wedge1].
\]
In the sequel, we will make use of these notations,
\begin{itemize}
\item  $B$ denotes the $d$-dimensional canonical process, i.e. $B_t(\omega)=\omega_t$, for any $\omega\in\Omega$.
\item $\mathcal{B}(\Omega)$ denotes the Borel $\sigma$-algebra of $\Omega$, similarly we have $\mathcal{B}(\Omega_{t})$ with  $\Omega_{t}:=\{ \omega_{\cdot \wedge t}:\omega \in \Omega \}$.
\item $L_{ip}(\Omega):=\{ \varphi(B_{t_{1}},\ldots,B_{t_{k}}):k\in \mathbb{N}
,t_{1},\ldots,t_{k}\in \lbrack0,\infty),\varphi \in C_{b.Lip}(\mathbb{R}
^{k\times d })\}$, where $C_{b.Lip}(\mathbb{R}^{k\times d})$ collects all
bounded Lipschitz functions on $\mathbb{R}^{k\times d}$, $L_{ip}(\Omega_{t})$ denotes all $\mathcal{B}(\Omega_{t})$-measurable elements in $L_{ip}
(\Omega)$.
\item $\mathbb{S}_{d}$ denote all symmetric matrices of size $d$.
\end{itemize}

For any given monotonic sublinear   continuous function $G:\mathbb{S}_d\rightarrow \mathbb{R}$, Peng \cite{P10} associated it with a nonlinear $G$-expectation $\mathbb{\hat{E}}[\cdot]$ using a nonlinear parabolic PDE,
which in turn makes the canonical processes $B$ a $d$-dimensional $G$-Brownian motion, ending up with the so-called $G$-expectation space $(\Omega,L_{ip}(\Omega),\mathbb{\hat{E}
}[\cdot],(\mathbb{\hat{E}}_{t}[\cdot])_{t\geq0})$. The readers are referred to \cite{P07,P08,P10} for detailed construction and so forth.

  For each
$p\geq1$, the completion of $L_{ip}(\Omega)$ under the norm
$||X||_{L_{G}^{p}}:=(\mathbb{\hat{E}}[|X|^{p}])^{1/p}$ is denoted by
$L_{G}^{p}(\Omega)$. Similarly, we can define $L_{G}^{p}(\Omega_t)$ for each $t>0$.
 The $G$-expectation $\mathbb{\hat{E}}[\cdot]$ and conditional $G$-expectation can be
extended continuously to the completion  $L_{G}^{p}(\Omega)$.
And the $G$-expectation can be regarded as a upper expectation.
\begin{theorem}[\cite{DHP2011,HP09}]
\label{the2.7}  There exists a weakly compact set
$\mathcal{P}$ of probability
measures on $(\Omega,\mathcal{B}(\Omega))$, such that
\[
\mathbb{\hat{E}}[\xi]=\sup_{P\in\mathcal{P}}E_{P}[\xi],\  \text{for
 all}\ \xi\in  {L}_{G}^{1}{(\Omega)}.
\]
\end{theorem}

For this $\mathcal{P}$, we define a capacity
\[
c(A):=\sup_{P\in\mathcal{P}}P(A),\ A\in\mathcal{B}(\Omega).
\]
A set $A\in\mathcal{B}(\Omega)$ is polar if $c(A)=0$.  A
property holds $``quasi$-$surely"$ (q.s.) if it holds except for a
polar set. In what follows,  we do not distinguish two random
variables between  $X$ and $Y$, if $X=Y$ q.s..

Now we state the nonlinear monotone convergence theorem, which is different from the linear case.
\begin{proposition}[\cite{DHP2011}]\label{downward convergence proposition} Suppose $X_n$, $n\geq 1$ and $X$ are $\mathcal{B}(\Omega)$-measurable.
If $\{X_n\}_{n=1}^\infty$ in ${L}_{G}^{1}(\Omega)$ satisfies that $X_n\downarrow X$, q.s.,  then
$\mathbb{\hat{E}}[X_n]\downarrow\mathbb{\hat{E}}[X].	$
\end{proposition}

Peng also introduced the stochastic integral with respect to $G$-Brownian motion, which led to a symmetric $G$-martingale.
Given a fixed contant $T>0$, the following spaces of stochastic processes will be useful,
\begin{itemize}
 \item $M_{G}^{0}(0,T):=\{\eta:\ \eta_{t}(\omega)=\sum\limits_{j=0}^{N-1}\xi_{j}(\omega)I_{[t_{j},t_{j+1})}(t),\  \xi_{i}\in L_{ip}(\Omega_{t_{i}})$ for some partition $t_{0}\leq t_1\leq \ldots\leq t_{N}$ of $[0,T]\}.$
 \item $M_{G}^{2}(0,T)$ is the completion of  $M_{G}^{0}(0,T)$ under norm $\Vert \eta \Vert_{M_{G}^{2}}=\{ \mathbb{\hat{E}}[\int_{0}^{T}|\eta_{s}|^{2}ds]\}^{1/2}.$
 \item $S_{G}^{0}(0,T)=\{h(t,B_{t_{1}\wedge t},\cdot \cdot \cdot,B_{t_{n}\wedge t}):t_{1},\ldots,t_{n}\in \lbrack0,T],h\in C_{b,Lip}(\mathbb{R}^{n+1})\}$.
 \item $S_{G}^{2}(0,T)$ is the completion of $S_{G}^{0}(0,T)$ under the norm $\Vert \eta\Vert_{S_{G}^{2}}=\{ \mathbb{\hat{E}}[\sup_{t\in\lbrack0,T]}|\eta_{t}|^{2}]\}^{\frac{1}{2}}$.
\end{itemize}

For each $1\leq i, j \leq d$,   we denote by $\langle B^i, B^j\rangle$   the mutual  variation process.
Then for two processes $ \eta\in M_{G}^{2}(0,T)$ and $ \xi\in M_{G}^{1}(0,T)$,
the $G$-It\^{o} integrals  $\int^{\cdot}_0\eta_sdB^i_s$ and $\int^{\cdot}_0\xi_sd\langle
B^i,B^j\rangle_s$  are well defined, see Peng \cite{P10}.
Moreover,   the corresponding $G$-It\^{o} formula were established.

\subsection{ $G$-BSDEs with Lipschitz assumptions}
From now on, we always assume that the function $G$  is  non-degenerate throughout our paper, i.e. there are two constants $0<\underline{\sigma}^{2}\leq\bar{\sigma}^{2}<\infty$ such that
$$ \frac{1}{2}\underline{\sigma}^2\text{tr}[A-A']\leq G(A)-G(A')\leq \frac{1}{2}\bar{\sigma}^2\text{tr}[A-A'], \text{ for } A\geq A'.$$
Then there exists a bounded and closed subset $\Gamma
\subset\mathbb{S}^+(d)$ such that
\begin{equation*}
G(A)=\frac{1}{2}\underset{Q\in\Gamma}{\sup}\mathrm{tr}[AQ],
\label{G-representation}
\end{equation*}
where $\mathbb{S}^+(d)$ denotes the space of all $d\times d$ symmetric positive-definite matrices.

Consider the following $G$-BSDEs (recall  that we use Einstein summation convention):
\begin{equation}\label{GBSDE}
 Y_t=\xi +\int_t^T f(s,Y_s,Z_s)ds+\int_{t}^{T}g^{ij}(s,Y_s,Z_s)d\langle B^i, B^j\rangle_{s}-\int_t^T  Z_s dB_s- (K_T-K_t),
 \end{equation}
in which the generators $$f(t,\omega,y,z), g^{ij}(t,\omega,y,z): [0,T]\times\Omega_T\times\mathbb{R}\times\mathbb{R}^d\rightarrow\mathbb{R}$$ and the terminal value $\xi$ fulfill these assumptions,
\begin{flalign*}
&\textbf{(A1)} \text{ For some }  \beta >2, \xi\in L_G^\beta(\Omega_T); \text{ for any }(y,z), f(\cdot, \cdot,y,z)\in M_G^\beta(0,T),g^{ij}(\cdot, \cdot,y,z)\in M_G^\beta(0,T).   &     \\
&\textbf{(A2)} \text{ There is a Lipschitz constant } L_0>0, \text{so that}    &
\end{flalign*}
$$ \vert f(t,y,z)-f(t,y',z')\vert +\vert g^{ij}(t,y,z)-g^{ij}(t,y',z')\vert  \leq L_0\vert y-y'\vert +L_0 \vert z-z'\vert.$$

For simplicity, we denote by $\mathfrak{S}_G^{2}(0,T)$ the collection of process $(Y,Z,K)$ such that $Y\in S_G^{2}(0,T)$, $Z\in M_G^{2}(0,T;\mathbb{R}^d)$, $K$ is a non-increasing $G$-martingale with $K_0=0$ and $K_T\in L_G^{2}(\Omega_T)$. Hu et al. \cite{HJPS2014,HJPS2014C} firstly obtained the existence and uniqueness result on Lipschitz $G$-BSDEs \eqref{GBSDE}, and the  comparison principle.
\begin{theorem}[\cite{HJPS2014}]\label{G1}
 Assume the conditions $\textbf{(A1)}$ and $\textbf{(A2)}$ hold.
   Then the equation $(\ref{GBSDE})$ admits a unique solution  $(Y,Z,K)\in \mathfrak{S}_G^{2}(0,T)$.
 \end{theorem}
\begin{theorem}[\cite{HJPS2014C}]\label{my2}
Assume $(\xi^{\nu},f^{\nu},g^{\nu}_{ij})$ satisfy assumption $\textbf{(A1)}$  for $\nu=1,2$.
Moreover,  one of them  satisfies assumption $\textbf{(A2)}$.
Suppose $(Y^{\nu},Z^{\nu},K^{\nu})$ is a $\mathfrak{S}_G^{2}(0,T)$-solution to the $G$-BSDE \eqref{GBSDE} with data $(\xi^{\nu},f^{\nu},g^{\nu}_{ij})$.
 If $\xi^2\leq \xi^1,$ $f^2\leq f^1$ and the matrix $(g^2_{ij})_{i,j=1}^d\leq (g^1_{ij})_{i,j=1}^d$, then we have $Y_t^2\leq Y_t^1$ for all $t\in [0,T].$
\end{theorem}

The     linear   $G$-BSDEs will be repeatedly used in our paper, so we sketch the idea on how to construct the solution. Consider linear $G$-BSDE of the form,
\begin{equation}\label{1LBSDE1}
 Y_t=\xi +\int_t^T [a_{s}Y_{s}+b_{s}Z_{s}+m_{s}]ds+\int_{t}^{T}[c_{s}^{ij}Y_{s}+d_{s}^{ij}Z_{s}+n^{ij}_{s}]d\langle B^i,B^j\rangle_{s}-\int_t^T  Z_s dB_s- (K_T-K_t),
 \end{equation}
  where $(a_{s})_{s\in\lbrack0,T]}$,
$(c^{ij}_{s})_{ s\in\lbrack0,T]}\in M_{G}^{2}(0,T)$, $(b_{s})_{s\in\lbrack0,T]}$,$(d^{ij}_{s})_{s\in\lbrack0,T]}\in M_{G}^{2}(0,T;\mathbb{R}^d)$ are bounded
processes and $\xi\in L_{G}^{2}(\Omega_{T})$,
$(m_{s})_{s\in\lbrack0,T]}$, $(n^{ij}_{s})_{s\in\lbrack0,T]}\in M_{G}^{2}(0,T)$.

To find the closed-form solution  to equation (\ref{1LBSDE1}), a standard method is to introduce a dual process. However for  the $G$-expectation case, unless the $G$-Brownian motion degenerates to the standard Wiener process,  the measures $ds$ and $d\langle B\rangle_{s}$ are mutually singular, therefore to cancel terms involving  $ds$ and $d\langle B^i,B^j\rangle_{s}$  is even harder.
To adapt the classical dual method, Hu et al. \cite{HJPS2014C} came up with a strategy of enlarging the original $G$-expectation space to
$\tilde{G}$-expectation space $(\tilde{\Omega},L_{\tilde{G}}^{1}
(\tilde{\Omega}),\mathbb{\hat{E}}^{\tilde{G}})$ with $\tilde{\Omega}=C_{0}([0,\infty),\mathbb{R}^{2d})$ and
\begin{align}\label{yw7}
\tilde{G}(A)=\frac{1}{2}\sup_{Q\in\Gamma}\mathrm{tr}\left[  A\left[
\begin{array}
[c]{cc}
Q & I_d\\
I_d & Q^{-1}
\end{array}
\right]  \right]  ,\ A\in\mathbb{S}_{2d}.
\end{align}
Let $(B_{t},\tilde{B}_{t})_{t\geq 0}$ be the canonical process in the extended space. Then
\begin{lemma}[\cite{HJPS2014C}]
\label{the5.2} In the extended $\tilde{G}$-expectation space, the solution of
the linear $G$-BSDE (\ref{1LBSDE1}) can be represented as
\begin{equation*}
Y_{t}=\mathbb{\hat{E}}_{t}^{\tilde{G}}[\tilde{\Gamma}_{T}^t\xi+\int_{t}^{T}
m_{s}\tilde{\Gamma}_s^tds+\int_{t}^{T}n_{s}^{ij}\tilde{\Gamma}^t_sd\langle B^i,B^j\rangle_{s}],
\end{equation*}
where $\{\tilde{\Gamma}^t_s\}_{s\in\lbrack t,T]}$ is the solution of the  following $\tilde{G}$-SDE:\begin{equation}
\tilde{\Gamma}^t_s=1+\int_{t}^{s}a_{r}\tilde{\Gamma}^t_rdr+\int_{t}^{s}c^{ij}_{r}\tilde{\Gamma}^t_rd\langle B^i, B^j\rangle
_{r}+\int_{t}^{s}d_{r}^{ij}\tilde{\Gamma}^t_rdB_{r}+\int_{t}^{s}b_{r}\tilde{\Gamma}^t_{r}d\tilde{B}_{r}.
\label{LSDE2}
\end{equation}
Moreover,
\begin{equation}\label{yw2}
\mathbb{\hat{E}}_{t}^{\tilde{G}}[\tilde{\Gamma}^t_TK_{T}-\int_{t}^{T}a_{s}K_{s}\tilde{\Gamma}^t_s
ds-\int_{t}^{T}c^{ij}_{s}K_{s}\tilde{\Gamma}^t_sd\langle
B^i, B^j\rangle_{s}]=K_{t}.
\end{equation}
\end{lemma}

\section{$G$-BSDEs with uniformly continuous generators}
In this section, we  shall investigate  the well-posedness problem of the subsequent $G$-BSDEs
\begin{equation}\label{GBSDE1}
 Y_t=\xi +\int_t^T f(s,Y_s,Z_s)ds+\int_t^T g^{ij}(s,Y_s,Z_s)d\langle B^i,B^j\rangle_s-\int_t^T  Z_s dB_s- (K_T-K_t),
 \end{equation}
 where the generators
\begin{equation*}
f(t,\omega,y,z),\  g^{ij}(t,\omega,y,z):[0,T]\times \Omega\times \mathbb{R}\times\mathbb{R}^d\rightarrow
\mathbb{R},
\end{equation*}
satisfy the following assumptions:
\begin{description}
\item[(H1)]  There exists a constant $ \beta >2 $ such that $f(\cdot, \cdot,y,z), g(\cdot,\cdot,y,z) \in M_G^\beta(0,T)$  for any $ y,z$.
\item[(H2)]   $f$ {and} $ g $ are Lipschitz continuous in $ y$, are of linear growth and  uniformly
 continuous in $ z$, i.e. there is a constant $ L$ and a continuous function $\phi$,  both independent of $(t,\omega)$,  such that
   $$ \vert f(t,\omega,y,z)-f(t,\omega,y',z')\vert+\vert g^{ij}(t,\omega,y,z)-g^{ij}(t,\omega,y',z')\vert\leq L\vert y-y'\vert+\phi(\vert z-z'\vert), $$
   where $\phi:\mathbb{R^+}\rightarrow \mathbb{R^+}$ is nondecreasing and sub-additive, with $\phi(0)=0$ as well as
   $\phi(z)\leq L(1+\vert z\vert).$
\item[(H3)]  $g^{ij}\equiv 0$ whenever $i\neq j$.
\end{description}

\begin{remark}
{\upshape
Note that assumption (H3) is necessary to construct a sequence of $G$-BSDEs monotonically converges to $Y$, see (i) of Lemma \ref{lemma2}.
}
\end{remark}

According to  Lemma 1 in Lepeltier and Martin \cite{LM} or Lemma 2 in Jia \cite{J},
there exists a sequence of  Lipschitz functions that nicely approximates $f$ and $g_{ij}$ respectively. Indeed,
for  any $(t,y,z),n\in\mathbb{N}$ and for every $\omega$, denote
\begin{align*}
\underline{\varphi}_n(t,y,z):=\inf\limits_{q\in\mathbb{Q}}\{\varphi(t,y,q)+n\vert z-q\vert\}-\varphi_0(t),\ \
\bar{\varphi}_n(t,y,z):=\sup\limits_{q\in\mathbb{Q}}\{\varphi(t,y,q)-n\vert z-q\vert\}-\varphi_0(t),
\end{align*}
where $\varphi=f, g^{ij}$ and $\varphi_0(t)=\varphi(t,0,0)$. Their main technical lemma can be summarized  as,
\begin{lemma}
\label{lemma1}
 Assume \textbf{(H1)-(H2)} hold. Then  for each $n>L$, the following properties hold
 \begin{description}
\item[(i)]  both $ \underline{\varphi}_n$ and $\bar{\varphi}_n$ are of linear growth, moreover, for all $(t,y,z)$,
    $$ -L(1+\vert y\vert+\vert z\vert)\leq \underline{\varphi}_n(t,y,z)\leq \varphi(t,y,z)-\varphi_0(t)\leq \bar{\varphi}_n(t,y,z) \leq  L(1+\vert y\vert+\vert z\vert);$$
\item[(ii)]  for all $(t,y,z)$,   $\underline{\varphi}_.(t,y,z)$ is  non-decreasing and $\bar{\varphi}_.(t,y,z)$ is non-increasing;
\item[(iii)]$\underline{\varphi}_n(t,y,\cdot)$ and $\bar{\varphi}_n(t,y,\cdot)$ are Lipschitz functions with constant $n$,
$\underline{\varphi}_n(t,\cdot,z)$ and $\bar{\varphi}_n(t,\cdot,z)$ are Lipschitz functions with Lipschitz constant $L$;
\item[(iv)] if $(y_n,z_n)\rightarrow (y,z),$  then $\underline{\varphi}_n(t,y_n,z_n)\rightarrow \varphi(t,y,z)-\varphi_0(t)$ and $\bar{\varphi}_n(t,y_n,z_n)\rightarrow \varphi(t,y,z)-\varphi_0(t)$;
\item[(v)] for all $(t,\omega,y,z)$,
$$0\leq \varphi(t,y,z)-\varphi_0(t)-\underline{\varphi}_n(t,y,z)\leq \phi(\frac{2L}{n-L}) , 0\leq \bar{\varphi}_n(t,y,z)+\varphi_0(t)-\varphi(t,y,z)\leq \phi(\frac{2L}{n-L}). $$
\end{description}
\end{lemma}

Based on the above  approximation results, we  construct two sequences of $G$-BSDEs  corresponding respectively to $(\underline{f}_n,\underline{g}_{n}^{ij})$ and $(\bar{f}_n,\bar{g}_{n}^{ij})$, i.e.,
\begin{align}\label{my1}
\begin{split}
\underline{Y}^n_t=&\xi +\int_t^T [\underline{f}_n(s,\underline{Y}^n_s, \underline{Z}^n_s)+f_0(s)]ds+\int_t^T[ \underline{g}^{ij}_n(s,\underline{Y}^n_s, \underline{Z}^n_s)+g^{ij}_0(s)]d\langle B^i,B^j\rangle_s -\int_t^T  \underline{Z}^n_s dB_s\\
& \ \ \ \ \ - (\underline{K}^n_T-\underline{K}^n_t),\\
 \bar{Y}^n_t=&\xi +\int_t^T [\bar{f}_n(s,\bar{Y}^n_s, \bar{Z}^n_s)+f_0(s)]ds+\int_t^T [\bar{g}_n^{ij}(s,\bar{Y}^n_s, \bar{Z}^n_s)+g^{ij}_0(s)]d\langle B^i,B^j\rangle_s-\int_t^T  \bar{Z}^n_s dB_s \\
& \ \ \ \ \ -(\bar{K}^n_T-\bar{K}^n_t).
\end{split}
\end{align}

We need an additional assumption to ensure the existence of $\underline{Y}^n$ and $ \bar{Y}^n$:
\begin{description}
\item[(H4)] For each $n$ and for any $(y,z)$, $\underline{\varphi}_n(t,y,z)$ and $\bar{\varphi}_n(t,y,z)$ all belong to  $M^{\beta}_G(0,T)$, with $\varphi=f,g^{ij}$.
\end{description}

\begin{remark}
{\upshape As can be easily seen, assumption \textbf{(H4)} is imposed mainly to keep all processes under investigation lying in   space $M^{\beta}_G(0,T)$.
This condition can be verified for lots of situations.  For instance, assume  \textbf{ (H1)} hold.
Suppose for $\varphi=f,g$ that $\varphi(\cdot, \cdot, y, z)$ is uniformly continuous in $(t, \omega)$ with the modulus of continuity
  independent of $(y, z)$,
\[
\vert \varphi(t,\omega,y,z)-\varphi(t',\omega',y,z)\vert\leq \phi(|t-t'|+\sup\limits_{s\in[0,t]}|\omega(s)-\omega'(s)|).
\]
Then it is straightforward  to observe that $\underline{\varphi}_n(\cdot,\cdot,y,z)$ and $\bar{\varphi}_n(\cdot,\cdot, y,z)$ are uniformly continuous in $(t,\omega)$.
Recalling the property (i) from  Lemma \ref{lemma1}, we have $\underline{\varphi}_n(\cdot,\cdot,y,z)$ and $\bar{\varphi}_n(\cdot,\cdot, y,z)$ are bounded  and then \[
\lim\limits_{N\rightarrow\infty}\mathbb{\hat{E}}[\int^T_0|\underline{\varphi}_n(t,y,z)|^{\beta}\mathbf{1}_{\{|\underline{\varphi}_n(t,y,z)|\geq N\}}dt]=0, \ \ \lim\limits_{N\rightarrow\infty}\mathbb{\hat{E}}[\int^T_0|\bar{\varphi}_n(t,y,z)|^{\beta}\mathbf{1}_{\{|\bar{\varphi}_n(t,y,z)|\geq N\}}dt]=0.
\]
Thus by Theorem 4.16 in \cite{HWZ2016}, we know $\varphi$ satisfies assumption \textbf{(H4)}.
}
\end{remark}

The following lemma is important in our future discussion.

\begin{lemma}\label{lemma2}
Let $\xi$ be in $L^{\beta}_G(\Omega_T)$ and  the assumptions \textbf{ (H1)-(H4)} hold. Then  the $G$-BSDE \eqref{my1} has a unique $\mathfrak{S}_G^{2}(0,T)$-solution. Moreover, we have
\begin{description}
\item[(i)]  for any $n,m\in\mathbb{N}$,   the comparisons $\underline{Y}^n\leq\underline{Y}^{n+1}\leq\bar{Y}^{m+1}\leq\bar{Y}^m$  hold;
\item[(ii)]  both $\underline{Y}^n$ and $\bar{Y}^n$ are uniformly bounded in $S_G^2(0,T)$;
\item[(iii)]  for each $n>L$,
the differences between $\underline{Y}^n$ and $\bar{Y}^n  $ can be uniformly controlled, that is,
$$\vert \underline{Y}^n_t-\bar{Y}^n_t\vert\leq C_G\phi(\frac{2L}{n-L}), \ \forall t\in[0,T],$$
where $C_G$ is a constant depending on $G,L$ and $T$.
\end{description}
\end{lemma}
\begin{proof}
The proof is built on the conclusions of lemma \ref{lemma1}. By assumption \textbf{(H3)}, we have $\underline{g}^{ij}_n=\bar{g}^{ij}_n=0$ whenever $i\neq j$.
Thus $(\underline{g}^{ij}_n-\underline{g}^{ij}_{n+1})_{i,j=1}^d$ is a nonnegative definite matrix.
Then from Theorem \ref{G1} and the comparison theorem \ref{my2}, it is trivial to verify  (i)   in  view of assertions (i)-(ii) from lemma \ref{lemma1}.

In order to prove (ii), setting $w(y,z)=L(1+|y|+|z|)$,  consider the following $G$-BSDEs
\begin{align*}
&U_t=\xi +\int_t^T [w( U_s, V_s)+f_0(s)] ds+\int_t^T [w( U_s, V_s)+g_0^{ii}(s)] d\langle B^i,B^i\rangle_s-\int_t^T  V_s dB_s- (R_T-R_t),\\
&U'_t=\xi +\int_t^T [-w( U'_s, V'_s)+f_0(s)] ds+\int_t^T [-w( U'_s, V'_s)+g_0^{ii}(s)] d\langle B^i,B^i\rangle_s-\int_t^T  V'_s dB_s- (R'_T-R'_t).
\end{align*}
It follows from Theorem \ref{G1} that the above $G$-BSDE admits a unique $\mathfrak{S}_G^{2}(0,T)$-solution $(U,V,R)$ and $(U',V',R')$, respectively.
Then by  (i) of lemma \ref{lemma1} and the comparison theorem \ref{my2},  it holds that for any $n\in\mathbb{N}$
\[
U'_t\leq \underline{Y}^n_t\leq\bar{Y}^n_t \leq U_t, \ \forall t\in[0,T], \]
which implies the desired result.

Finally, we proceed to verify the third assertion (iii).  Without loss of generality, assume  that  $d=1$.
Set $(\hat{Y},\hat{Z})=(\bar{Y}^n-\underline{Y}^n,\bar{Z}^n-\underline{Z}^n)$.
Then for each $t\in[0,T]$, we have
 \begin{equation}\label{eq5}
\hat{Y}_t+\underline{K}^n_t=\underline{K}^n_T+\int_t^T \hat{f}_sds+\int^T_t\hat{g}_sd\langle B\rangle_s-\int_t^T  \hat{Z}_s  dB_s- (\bar{K}^n_T-\bar{K}^n_t),
 \end{equation}
where $ \hat{\varphi}_s=\bar{\varphi}_n(s,\bar{Y}^n_s, \bar{Z}^n_s)-\underline{\varphi}_n(s,\underline{Y}^n_s, \underline{Z}^n_s)$ for $\varphi=f,g.$

By Lemma 3.5 in \cite{HW}, for each $\varepsilon>0$, there exist four bounded processes $a^{\varepsilon}, b^{\varepsilon}$, $c^{\varepsilon}$, $d^{\varepsilon}\in M^2_G(0,T)$
such that for all $s\in\lbrack 0,T\rbrack,$
\[
\hat{f}_{s}=a_{s}^{\varepsilon}\hat{Y}_{s}+b_{s}^{\varepsilon}\hat{Z}
_{s}+m_s-m_{s}^{\varepsilon},\ \hat{g}_{s}=c_{s}^{\varepsilon}\hat{Y}
_{s}+d_{s}^{\varepsilon}\hat{Z}_{s}+n_s-n_{s}^{\varepsilon},
\]
and $|a_{s}^{\varepsilon}|\leq L$, $|c_{s}^{\varepsilon}|\leq L$,
$|b_{s}^{\varepsilon}|\leq n$, $|d_{s}^{\varepsilon}|\leq n$, $|m_{s}^{\varepsilon}|\leq2(L+n)\varepsilon$, $|n_{s}^{\varepsilon}
|\leq2(L+n)\varepsilon$, $m_{s}=\bar{f}_n(s,\underline{Y}^n_s, \underline{Z}^n_s)-\underline{f}_n(s,\underline{Y}^n_s, \underline{Z}^n_s)$,
$n_{s}=\bar{g}_n(s,\underline{Y}^n_s, \underline{Z}^n_s)-\underline{g}_n(s,\underline{Y}^n_s, \underline{Z}^n_s)$.

In order to estimate the solution to  the above linearized equation (\ref{eq5}), as in \cite{HJPS2014C}, we
 shift from the underlying $G$-expectation space to an auxiliary extended
$\tilde{G}$-expectation space $(\tilde{\Omega},L_{\tilde{G}}^{1}
(\tilde{\Omega}),\mathbb{\hat{E}}^{\tilde{G}})$ with $\tilde{\Omega}=C_{0}([0,\infty),\mathbb{R}^{2})$, where $\tilde{G}$ is given by equation \eqref{yw7}, within which,
  $(B_{t},\tilde{B}_{t})_{t\geq 0}$  denotes the corresponding canonical process.

Applying Lemma \ref{the5.2}  yields that
\begin{align*}
 \hat{Y}_{t}+\underline{K}^n_t
&  =\mathbb{\hat{E}}_{t}^{\tilde{G}}[\tilde{\Gamma}_{T}
^{t,\varepsilon}\underline{K}_{T}^{n}+\int_{t}^{T}(m_s+2G(n_s)-m_{s}^{\varepsilon
}-a_{s}^{\varepsilon}\underline{K}^n_s)\tilde{\Gamma}_{s}^{t,\varepsilon}ds\\
&  -\int_{t}^{T}(n_{s}^{\varepsilon}+c_{s}^{\varepsilon}\underline{K}^n_s
)\tilde{\Gamma}_{s}^{t,\varepsilon}d\langle B\rangle_{s}+\int_{t}^{T}n_s\tilde{\Gamma}_{s}^{t,\varepsilon}d\langle B\rangle_s-\int_{t}^{T}2G(n_s)\tilde{\Gamma}_{s}^{t,\varepsilon}ds],
\end{align*}
where $\{\tilde{\Gamma}^{t,\varepsilon}_{s}\}_{s\in\lbrack t,T\rbrack}$ is given by, c.f. equation \eqref{LSDE2}. From $G$-It\^{o}'s formula, we conclude  that
\begin{equation*}
\tilde{\Gamma}^{t,\varepsilon}_{s}=\exp(\int_{t}^{s}(a^{\varepsilon}_{r}-b^{\varepsilon}_{r}d^{\varepsilon}_{r})dr+\int_{t}^{s}c^{\varepsilon}_{r}d\langle
B\rangle_{r})\mathcal{E}_{s}^{B}\mathcal{E}_{s}^{\tilde{B}}.
\end{equation*}
Here $\mathcal{E}_{s}^{B}=\exp(\int^s_td^{\varepsilon}_{r}dB_r-\frac{1}{2}\int^s_t|d^{\varepsilon}_{r}|^2d\langle B\rangle_r)$ and
$\mathcal{E}_{s}^{\tilde{B}}=\exp(\int^s_tb^{\varepsilon}_{r}d\tilde{B}_r-\frac{1}{2}\int^s_t|b^{\varepsilon}_{r}|^2d\langle \tilde{B}\rangle_r)$.
Therefore using equations \eqref{eq5} and $G$-It\^{o}'s formula we get that
\begin{align}
 \hat{Y}_{t}+\underline{K}^n_t
  \leq \mathbb{\hat{E}}_{t}^{\tilde{G}}
[\int_{t}^{T}(m_s+2G(n_s))\tilde{\Gamma}_{s}^{t,\varepsilon}ds-\int_{t}^{T}m_{s}^{\varepsilon}\tilde{\Gamma}_{s}^{t,\varepsilon}ds
-\int_{t}^{T}n_{s}^{\varepsilon
}\tilde{\Gamma}_{s}^{t,\varepsilon}d\langle B\rangle_{s}]+\underline{K}^n_t,\ \ q.s..\label{myq4}
\end{align}
By (v) of Lemma \ref{lemma1}, we get
\[
0\leq m_s+2G(n_s)\leq 2(1+\bar{\sigma}^2)\phi(\frac{2L}{n-L}).
\]
Note that for each $s\geq t$, $
\tilde{\Gamma}^{t,\varepsilon}_s\leq \exp(L(1+\bar{\sigma}^2)(s-t)){\Gamma}^{t,\varepsilon}_s,
$
where $
{\Gamma}^{t,\varepsilon}_{s}=1+\int_{t}^{s}d^{\varepsilon}_{r}{\Gamma}^{t,\varepsilon}_{r}dB_{r}+\int_{t}
^{s}b^{\varepsilon}_{r}{\Gamma}^{t,\varepsilon}_{r}d\tilde{B}_{r}.
$
Then by equation \eqref{myq4},
we derive that
\begin{align*}
\hat{Y}_{t}\leq &
[2(1+\bar{\sigma}^2)\phi(\frac{2L}{n-L})+2(L+n)(1+\bar{\sigma}^2)\varepsilon]\mathbb{\hat{E}}^{\tilde{G}}_t[\int_{t}^{T}\exp(L(1+\bar{\sigma}^2)(s-t)){\Gamma}^{t,\varepsilon}_{s}ds]\\
\leq &\frac{\exp(L(1+\bar{\sigma}^2)(T-t))}{L}[2\phi(\frac{2L}{n-L})+2(L+n)\varepsilon].
\end{align*}
Sending $\varepsilon\rightarrow 0$,  we have
\[
\hat{Y}_{t} \leq \frac{2\exp(L(1+\bar{\sigma}^2)(T-t))}{L}\phi(\frac{2L}{n-L}),
\]
which completes the proof.
\end{proof}
\begin{remark}{\upshape
Note that from (i) and (ii) of Lemma \ref{lemma2}, in general  we cannot conclude that $\underline{Y}^n$ (or $\bar{Y}^n$) is a Cauchy sequence in $M^{2}_G(0,T)$  according to Proposition \ref{downward convergence proposition},
which is different from the classical case.
}
\end{remark}

Now we are ready to state the main result of this section.

\begin{theorem}\label{main}
Given assumptions \textbf{(H1)-(H4)} and $\xi\in L^{\beta}_G(\Omega_T)$, the $G$-BSDE (\ref{GBSDE1}) admits a unique  solution  $(Y,Z,K)\in \mathfrak{S}_G^{2}(0,T).$
\end{theorem}

\begin{proof}
We shall deal with the existence and uniqueness of solution to $G$-BSDE (\ref{GBSDE1}) separately.  For the  uniqueness, suppose that  both of $(Y^i,Z^i,K^i)$, $i=1,2$ are $\mathfrak{S}_G^{2}(0,T)$-solution to $G$-BSDE (\ref{GBSDE1}), by comparison theorem \ref{my2}, we obtain that for each $n$
\[
\underline{Y}^n_t\leq Y^i_t\leq \bar{Y}^n_t, \ \forall t\in [0,T],
\]
which, together with Lemma \ref{lemma2}, implies
\[
|Y^1_t-Y^2_t|\leq |\bar{Y}^n_t-\underline{Y}^n_t| \leq C_G\phi(\frac{2L}{n-L}), \ \forall n>L.
\]
Note that $Y^i_t$ is a continuous process.
Sending $n\rightarrow\infty$, we deduce that $Y^1=Y^2$ q.s.. Then applying $G$-It\^{o}'s formula upon $\vert Y^1_s-Y^2_s\vert^2$ on $[0,T]$,  we have
$Z^1=Z^2$ and then $K^1=K^2$, which shows that $G$-BSDE (\ref{GBSDE1}) allows for at most one $\mathfrak{S}_G^{2}(0,T)$-solution.

The rest of the proof is devoted  to studying the existence, which will be divided into  three steps. Without loss of generality, we assume $d=1$ and  $g\equiv 0$.

{\it 1 The uniform estimates.}
Let $C(\alpha)$ denote a    constant depending on  parameter $\alpha$  that may change from line to line.
From (ii) of Lemma \ref{lemma2}, we have for all $n$ \[
\|\bar{Y}^n \|_{S^2_G}\leq C(L, \bar{\sigma},\underline{\sigma}, T).\]

Calculating by It\^{o}'s formula upon $\vert \bar{Y}^n \vert^2$, we have for any $t\in[0,T]$,
\begin{align}\label{my3}
\vert \bar{Y}^n_t\vert^2+\int_t^T \vert  \bar{Z}^n_s\vert^2 d\langle B\rangle_s=\vert \xi\vert^2+2\int_t^T\bar{Y}^n_s(\bar{f}_n(s,\bar{Y}^n_s, \bar{Z}^n_s)+f_0(s))ds-2\int_t^T\bar{Y}^n_s \bar{Z}^n_sdB_s-\int_t^T2\bar{Y}^n_sd\bar{K}^n_s.
\end{align}
Since \[
 |\bar{f}^n(t,y,z)|\leq L(1+|y|+|z|),
 \]
we get that
\[
2\bar{Y}^n_s\bar{f}_n(s,\bar{Y}^n_s, \bar{Z}^n_s)\leq 2L(|\bar{Y}^n_s|+|\bar{Y}^n_s|^2)+\frac{4L^2}{\underline{\sigma}^2}|\bar{Y}^n_s|^2+\frac{\underline{\sigma}^2}{4}|\bar{Z}^n_s|^2.
\]
Using BDG inequality and H\"{o}lder's inequality, we derive that
\begin{align*}
 \mathbb{\hat{E}}[\vert \int_0^T \bar{Y}^n_s \bar{Z}^n_sdB_s\vert]\leq C(\bar{\sigma})\mathbb{\hat{E}}[\vert \int_0^T |\bar{Y}^n_s \bar{Z}^n_s|^2ds\vert^{\frac{1}{2}}] &\leq C(\bar{\sigma})\|\bar{Y}^n_s\|_{S^2_G}\|\bar{Z}^n_s\|_{M^2_G}\\
 &\leq C(\bar{\sigma},\underline{\sigma})\|\bar{Y}^n_s\|_{S^2_G}^2+\frac{\underline{\sigma}^2}{8}\|\bar{Z}^n_s\|_{M^2_G}.
\end{align*}
Thus, in view of equation \eqref{my3} we have
\begin{align}\label{my4}
 \mathbb{\hat{E}}[\int_0^T \vert  \bar{Z}^n_s\vert^2 d\langle B\rangle_s]\leq  C(L, \bar{\sigma},\underline{\sigma}, T)+\frac{\underline{\sigma}^2}{2}\|\bar{Z}^n_s\|_{M^2_G}+
 2\mathbb{\hat{E}}[\sup\limits_{s\in[0,T]}|\bar{Y}^n_s||\bar{K}^n_T|].
\end{align}
Recalling that \[
\bar{K}^n_T=\xi-\bar{Y}^n_0+\int_0^T [\bar{f}_n(s,\bar{Y}^n_s, \bar{Z}^n_s)+f_0(s)]ds-\int_0^T  \bar{Z}^n_s dB_s.
\]
By a similar analysis as above,  we obtain
\[
\mathbb{\hat{E}}[\sup\limits_{s\in[0,T]}|\bar{Y}^n_s||\bar{K}^n_T|]\leq  C(L, \bar{\sigma},\underline{\sigma}, T)+\frac{\underline{\sigma}^2}{8}\|\bar{Z}^n_s\|_{M^2_G},
\]
putting together   equation \eqref{my4} with the fact that $\underline{\sigma}^2 \Vert \bar{Z}^n\Vert^2_{M_{G}^2(0,T)} \leq \mathbb{\hat{E}}[\int_0^T \vert  \bar{Z}^n_s\vert^2 d\langle B\rangle_s]$
indicates that
$$ \Vert \bar{Z}^n\Vert_{M_{G}^2(0,T)}+\|\bar{K}^n_T\|_{L^2_G} \leq  C(L, \bar{\sigma},\underline{\sigma}, T), \ \forall n\in\mathbb{N}.$$

{\it 2 The convergence.} From assertions  (i) and (iii) of Lemma \ref{lemma2}, we get that for each $n,m>L$
\[
\|\bar{Y}^n-\bar{Y}^m\|_{S_G^2}\leq \|\bar{Y}^{n\wedge m}-\underline{Y}^{n\wedge m}\|_{S_G^2}\leq  C_G\phi(\frac{2L}{n\wedge m -L})
\]
from which we conclude that $\{ \bar{Y}^n\}_{n\in\mathbb{N}}$  is a Cauchy sequence in $S_G^2(0,T)$.
Then there is a process $Y\in S_G^2(0,T)$ such that $\bar{Y}^n$ converges to $Y$ in $S_G^2(0,T)$.

We continue  to show the convergence of $\bar{Z}^n$ in $ M_G^2(0,T)$.
For each $n,m>L$, applying It\^{o}'s formula to $|\bar{Y}^n-\bar{Y}^m|^2$ yields that
\begin{align*}
&\underline{\sigma}^2\mathbb{\hat{E}}[\int_0^T \vert\bar{Z}^n_s-\bar{Z}^m_s\vert^2 ds]\leq \mathbb{\hat{E}}[\int_0^T \vert\bar{Z}^n_s-\bar{Z}^m_s\vert^2 d\langle B\rangle_s]\\
&\leq 2\mathbb{\hat{E}}[\int_0^T(\bar{Y}^n_s-\bar{Y}^m_s)(\bar{f}_n(s,\bar{Y}^n_s, \bar{Z}^n_s)-\bar{f}_m(s,\bar{Y}^m_s, \bar{Z}^m_s))ds
-\int_0^T(\bar{Y}^n_s-\bar{Y}^m_s)d(\bar{K}^n_s-\bar{K}^m_s)]\\
&\leq  2\mathbb{\hat{E}}[\sup\limits_{s\in[0,T]}|\bar{Y}^n_s-\bar{Y}^m_s|\{L\int_0^T(2+|\bar{Y}^n_s|+|\bar{Y}^m_s|+|\bar{Z}^n_s|+|\bar{Z}^n_s|)ds + |\bar{K}^n_T|+|\bar{K}^m_T|\}]\\
&\leq C(L, \bar{\sigma},\underline{\sigma}, T)\|\bar{Y}^n-\bar{Y}^m\|_{S^2_G},
\end{align*}
where we have used the estimates of {\it step 1} and H\"{o}lder's inequality in the last inequality.
Consequently, we can find some process $Z\in M^2_G(0,T)$ so that  $\bar{Z}^n$ converges to $Z$ in $M_G^2(0,T)$.

Denote
\[
{K}_t:=Y_t-Y_0+\int_0^tf(s,{Y}_s,{Z}_s)ds-\int_0^t{Z}_s dB_s,
\]
we claim that
\begin{align}\label{my5}
\lim_{n\rightarrow\infty}\mathbb{\hat{E}}[\int_{0}^{T}\vert \bar{f}_n(s,\bar{Y}^n_s, \bar{Z}^n_s)+f_0(s)-f(s,{Y}_s,{Z}_s) \vert^{2}ds]=0,
\end{align}
whose proof will be given in {\it step 3}.
Thus it is easy to check that for each $t\in[0,T]$
\[
\lim_{n\rightarrow\infty}\mathbb{\hat{E}}[|{K}_t-\bar{K}^n_t|^2]=0,
\]
which implies that $K$ is a non-increasing $G$-martingale and then $(Y,Z,K)\in\mathfrak{S}_G^{2}(0,T) $ is the solution to $G$-BSDE \eqref{GBSDE1}.

{\it 3 The proof of equation \eqref{my5}.} For each $n>L$, applying lemma \ref{lemma1}, we get that
\begin{align*}
  &\mathbb{\hat{E}}[\int_{0}^{T}\vert \bar{f}_n(s,\bar{Y}^n_s, \bar{Z}^n_s)+f_0(s)-f(s,{Y}_s,{Z}_s) \vert^{2}ds] \\
 & \leq  3\mathbb{\hat{E}}[\int_{0}^{T}\vert \bar{f}_n(s,\bar{Y}^n_s, \bar{Z}^n_s)+f_0(s)-f(s,\bar{Y}^n_s,\bar{Z}^n_s) \vert^{\alpha}ds +\int_{0}^{T}\vert f(s,\bar{Y}^n_s,\bar{Z}^n_s)-f(s,Y_s,\bar{Z}^n_s)\vert^2 ds\\
 &\ \ \ \ +\int_{0}^{T}\vert f(s,{Y}_s,\bar{Z}^n_s)-f(s,Y_s,{Z}_s)\vert^2 ds]
 \\
 & \leq 3 T \phi(\frac{2L}{n-L})+3T L^2 \|\bar{Y}^n-Y\|_{S^2_G}^2+3\mathbb{\hat{E}}[\int_{0}^{T}\vert f(s,{Y}_s,\bar{Z}^n_s)-f(s,Y_s,{Z}_s)\vert^2 ds].
\end{align*}
By the uniform continuity of $f$ in $z$,  for any fixed $\varepsilon >0$, there exists a $\delta$, so that $ \vert f(\cdot,\cdot,x)-f(\cdot,\cdot,y)\vert <\varepsilon$ whenever  $\vert z-q\vert \leq\delta$.
Then for each $N>0$, we obtain that
\begin{align*}
  &\mathbb{\hat{E}}[\int_{0}^{T}\vert f(s,{Y}_s,\bar{Z}^n_s)-f(s,Y_s,{Z}_s)\vert^2 ds]\\ &\leq
  2\mathbb{\hat{E}}[\int_{0}^{T}\vert f(s,{Y}_s,\bar{Z}^n_s)-f(s,Y_s,{Z}_s)\vert^2 \mathbf{1}_{\vert \bar{Z}^n_s-Z_s\vert\leq \delta} ds]
  +2\mathbb{\hat{E}}[\int_{0}^{T}\vert f(s,{Y}_s,\bar{Z}^n_s)-f(s,Y_s,{Z}_s)\vert^2 \mathbf{1}_{\vert \bar{Z}^n_s-Z_s\vert > \delta} ds] \\
   &\leq  2 T \varepsilon^2 + 2\mathbb{\hat{E}}[\int_{0}^{T}(|f_0(s)|+2L\vert Y_s\vert+L\vert \bar{Z}^n_s\vert +L\vert Z_s\vert )^2\mathbf{1}_{\vert \bar{Z}^n_s-Z_s\vert > \delta} ds],
\end{align*}
Since  $\bar{Z}^n$ converges to $Z$ in $M_G^2(0,T)$, it is easy to check that
$\mathbb{\hat{E}}[\int^T_0\mathbf{1}_{\vert \bar{Z}^n_s-Z_s\vert > \delta}]$ is vanishing as $n\rightarrow\infty$.
 Note that $|f_0(s)|+2L\vert Y_s\vert+L\vert \bar{Z}^n_s\vert +L\vert Z_s\vert\in M^2_G(0,T) $. Thus with the help of Theorem 4.7 in \cite{HWZ2016}, we get that
 \[
 \lim_{n\rightarrow\infty}\mathbb{\hat{E}}[\int_{0}^{T}(|f_0(s)|+2L\vert Y_s\vert+L\vert \bar{Z}^n_s\vert +L\vert Z_s\vert )^2\mathbf{1}_{\vert \bar{Z}^n_s-Z_s\vert > \delta} ds]=0.
 \]

 Consequently, putting together the above two inequalities we deduce that
 \[
 \limsup_{n\rightarrow\infty}\mathbb{\hat{E}}[\int_{0}^{T}\vert \bar{f}_n(s,\bar{Y}^n_s, \bar{Z}^n_s)+f_0(s)-f(s,{Y}_s,{Z}_s) \vert^{2}ds] \leq 2 T \varepsilon^2.
 \]
Letting $\varepsilon\rightarrow 0$, we get the desired result.
\end{proof}

\begin{example}{\upshape
For a $1$-dimensional $G$-Brownian motion $B$ with $\underline{\sigma}^2:= -\mathbb{\hat{E}}[-|B_1|^2]$, consider the following $G$-BSDE:
\[
Y_t=\frac{1}{6}|B_T|^6-\frac{5}{2}\underline{\sigma}^2\int^T_t|Z_s|^{\frac{4}{5}}ds-\int^T_t Z_sdB_s-(K_T-K_t).
\]
Note that $f(z)=-\frac{5}{2}\underline{\sigma}^2|z|^{\frac{4}{5}}$ is a uniformly continuous function.
Then by $G$-It\^{o}'s formula and Theorem \ref{main}, it is easy to check that $(\frac{1}{6}|B_t|^6,(B_t)^5,\frac{5}{2}\underline{\sigma}^2\int^t_0 |B_s|^4ds- \frac{5}{2}\int^t_0 |B_s|^4d\langle B\rangle_s)$
is the unique  $\mathfrak{S}_G^{2}(0,T)$-solution.
}
\end{example}

\begin{theorem}[Comparison Theorem]\label{HM7}
Suppose $\xi^\nu\in L^{\beta}_G(\Omega_T)$,  $\nu=1,2$ and $f^\nu,g^{\nu,ij}$ satisfy assumption (H1)-(H4).
Let $(Y^\nu,Z^\nu,K^\nu)$ be the $\mathfrak{S}_G^{2}(0,T)$-solution of $G$-BSDE \eqref{GBSDE1} with data  $(\xi^\nu,f^\nu,g^{\nu,ij})$.
If $\xi^1 \leq \xi^2$, $f^1(t,y,z)\leq f^2(t,y,z)$ and $g^{1,ij}(t,y,z)\leq g^{2,ij}(t,y,z)$ for any $(t,\omega,y,z)$,  then
$Y^1_t\leq Y^2_t$ for each $t$.
\end{theorem}
\begin{proof}
For each $n\in\mathbb{N}$, let $(\bar{Y}^{2,n},\bar{Z}^{2,n},\bar{K}^{2,n} )$ be the  $\mathfrak{S}_G^{2}(0,T)$-solution of $G$-BSDE \eqref{my1} corresponding to $(f^2,g^{2,ij})$.
It is obvious that $ \bar{f}^2_n(t,y,z)+f^2_0(t)\geq f^1(t,y,z)$ and $\bar{g}_n^{2,ii}(t,y,z)+g^{2,ii}_0(t)\geq g^{1,ii}(t,y,z)$. Note that $g^{2,ij}=g^{1,ij}=\bar{g}_n^{2,ij}=0$ whenever $i\neq j$.
Then using Theorem \ref{my2}, we get that ${Y}^1_t\leq \bar{Y}^{2,n}_t$ for each $t$.
Note that $\bar{Y}^{2,n}$ converges to $Y^2$ in $S^2_G(0,T)$.
Sending $n\rightarrow\infty$, we derive that $Y^1_t\leq Y^2_t$. The proof is complete.
\end{proof}

\section{Nonlinear  Feynman-Kac formula}
In this section, we shall utilize  Theorem  \ref{main} to establish a nonlinear  Feynman-Kac formula that slightly generalizes the corresponding result of \cite{HJPS2014C,P10}.
Retaining the notations in previous sections,  for each $(t,x)\in[0,T]\times\mathbb{R}^m$, let's consider    $G$-BSDE
\begin{align}
  dY_s^{t,x}=&-f(s,X_s^{t,x},Y_s^{t,x},Z_s^{t,x})ds-g^{i}(s,X_s^{t,x},Y_s^{t,x},Z_s^{t,x})d\langle B^i,B^i\rangle_s+Z_s^{t,x}dB_s+dK_s^{t,x}, \ s\in[t,T] \notag\\ Y_T^{t,x}=&\Phi(X_T^{t,x}),\label{fbsde2}
\end{align}
where $X^{t,x}$ is defined through  a forward $G$-SDE on the interval $[t,T]$
\begin{equation}\label{fbsde1}
  dX_s^{t,x}=b(s,X_s^{t,x})ds+h^{ij}(s,X_s^{t,x})d\langle B^i,B^j\rangle_s+\sigma(s,X_s^{t,x})dB_s, \ \ X_t^{t,x}=x.
\end{equation}

In the sequel, we use these   running assumptions abbreviated as \textbf{(H5)} :
\begin{itemize}
  \item[(i)] $b,h^{ij}=h^{ji}:[0,T]\times\mathbb{R}^m\rightarrow \mathbb{R}^m; \sigma:[0,T]\times\mathbb{R}^{m}\rightarrow \mathbb{R}^{m\times d}; f,g^{i}:[0,T]\times\mathbb{R}^m\times\mathbb{R}\times\mathbb{R}^d\rightarrow\mathbb{R};\Phi:\mathbb{R}^m\rightarrow \mathbb{R}, $ are all deterministic continuous functions.
  \item[(ii)]There exist two positive integers $q, L$ and a  modulus of continuity $\phi$ such that
  \begin{align*}
     &\vert b(t,x)-b(t,x')\vert+\sum_{i,j=1}^{d}\vert h^{ij}(t,x)-h^{ij}(t,x')\vert
    +\vert \sigma(t,x)-\sigma(t,x')\vert \leq L\vert x-x'\vert,\\
    &\vert\Phi(x)-\Phi(x')\vert +\vert f(t,x,y,z)-f(t,x',y',z')\vert+\sum_{i=1}^{d}\vert g^{i}(t,x,y,z)-g^{i}(t,x',y',z')\vert\\
   &\leq L(1+\vert x\vert^q+\vert x'\vert^q)\vert x-x'\vert+L\vert y-y'\vert+ \phi(\vert z-z'\vert).
  \end{align*}
\end{itemize}

To link the above $G$-BSDE system with  PDE, we need several estimates
 from \cite{HJPS2014C,P10},
\begin{lemma}\label{estimateX}
  Assuming  \textbf{(H5)}, for any $\delta \in[0,T-t]$, there exists a constant $C$ depending on $L',G,p,n,T$ such that
  \begin{align*}
    & \mathbb{\hat{E}}_t[\vert X_{t+\delta}^{t,x}\vert^p]\leq C(1+\vert x\vert^p),\\
&\mathbb{\hat{E}}_t[\vert X_{t+\delta}^{t,x}-X_{t+\delta}^{t,x'} \vert^p]\leq C\vert x-x'\vert^p,\\
&\mathbb{\hat{E}}_t\Big[\sup_{s\in[t, t+\delta]}\vert X_s^{t,x}-x\vert^p\Big]\leq C(1+\vert x\vert^p)\delta^{p/2}.
  \end{align*}
\end{lemma}

\begin{theorem}
Suppose  \textbf{(H5)} hold. Then $G$-BSDE \eqref{fbsde2} has a unique solution triplet  $(Y^{t,x},Z^{t,x},K^{t,x})\in \mathfrak{S}_G^{2}(t,T).$
\end{theorem}
\begin{proof}
By Lemma \ref{estimateX} and assumption \textbf{(H5)}, for each $p\geq 1,$
it is easy to get that $\Phi(X_T^{t,x})\in L^p_G(\Omega_T)$. The facts $f(.,X_.^{t,x},y,z)$, $g^{ij}(.,X_.^{t,x},y,z)\in M^p_G(t,T)$  follow from Theorem 4.16 in \cite{HWZ2016}. Therefore it suffices to verify conditions \textbf{(H4)}, before applying Theorem \ref{main} to complete the proof.

For  any $(t,x,y,z)\in [0,T]\times\mathbb{R}^m\times\mathbb{R}\times\mathbb{R}^d$ and $n\in\mathbb{N}$,
set
\begin{align*}
\underline{\varphi}_n(t,x,y,z)& =\inf\limits_{q\in\mathbb{Q}}\{\varphi(t,x,y,q)+n\vert z-q\vert\}-\varphi_0(t,x),\\
\bar{\varphi}_n(t,x,y,z)& =\sup\limits_{q\in\mathbb{Q}}\{\varphi(t,x,y,q)-n\vert z-q\vert\}-\varphi_0(t,x),
\end{align*}
for $\varphi=f, g^{i}$ and $\varphi_0(t,x)=\varphi(t,x,0,0)$.
By property (ii) of \textbf{(H5)}, we derive that
\[
|\underline{\varphi}_n(t,x,y,z)-\underline{\varphi}_n(t,x^{\prime},y,z)|+|\bar{\varphi}_n(t,x,y,z)-\bar{\varphi}_n(t,x^{\prime},y,z)|\leq 4L(1+\vert x\vert^q+\vert x'\vert^q)\vert x-x'\vert,
\]
which, together with Theorem 4.16 in \cite{HWZ2016} and Lemma \ref{estimateX},  implies that both
$\underline{\varphi}_n(.,X_.^{t,x},y,z)$ and $\bar{\varphi}_n(.,X_.^{t,x},y,z)$ belong to  $M^p_G(t,T)$ for each $p\geq 1.$
\end{proof}

Using the same notations appearing  in the above argument,
for each $(t,x)\in [0,T]\times\mathbb{R}^m$ and $n\in \mathbb{N}$, we consider a sequence of approximating $G$-BSDEs corresponding respectively to generators $(\underline{f}_n,\underline{g}_{n}^{ij})$ and $(\bar{f}_n,\bar{g}_{n}^{ij})$ on $[t,T]$,
\begin{align}
\begin{split}
\underline{Y}^{n,t,x}_s=&\Phi(X^{t,x}_T) +
\int_s^T [\underline{f}_n(r,X^{t,x}_r,\underline{Y}^{n,t,x}_r, \underline{Z}^{n,t,x}_r)+f_0(r,X^{t,x}_r)]dr-\int_s^T  \underline{Z}^{n,t,x}_r dB_r
 \\
& \ \ \ \ \ +\int_s^T[ \underline{g}^{i}_n(r,X^{t,x}_r,\underline{Y}^{n,t,x}_r, \underline{Z}^{n,t,x}_r)+g^{i}_0(r,X^{t,x}_r)]d\langle B^i,B^i\rangle_r - (\underline{K}^{n,t,x}_T-\underline{K}^{n,t,x}_s),\\
 \bar{Y}^{n,t,x}_s=&\Phi(X^{t,x}_T) +\int_s^T [\bar{f}_n(r,X^{t,x}_r,\bar{Y}^{n,t,x}_r, \bar{Z}^{n,t,x}_r)+f_0(r,X^{t,x}_r)]dr -\int_t^T  \bar{Z}^{n,t,x}_s dB_s
 \\
& \ \ \ \ \ +\int_s^T [\bar{g}_n^{i}(r,X^{t,x}_r,\bar{Y}^{n,t,x}_r, \bar{Z}^{n,t,x}_r)+g^{i}_0(r,X^{t,x}_r)]d\langle B^i,B^i\rangle_r -(\bar{K}^{n,t,x}_T-\bar{K}^{n,t,x}_s).
\end{split}
\end{align}
If we denote
$$ \underline{u}^n(t,x):=\underline{Y}^{n,t,x}_t, \ \ \bar{u}^n(t,x):=\bar{Y}^{n,t,x}_t, \ \ (t,x)\in [0,T]\times\mathbb{R}^m.$$
By Proposition 4.2 in \cite{HJPS2014C}, both $\underline{u}^n$ and $\bar{u}^n$ are continuous functions. Similarly we can define
$$ u(t,x):=Y^{t,x}_t, \ \ (t,x)\in [0,T]\times\mathbb{R}^m.$$
Clearly $u(t,x)$ is a well-defined deterministic function from the above theorem. And some regularity can be derived from that of $\underline{u}^n,\bar{u}^n$,
indeed we have
\begin{lemma} \label{my8}
Given assumption \textbf{(H5)}, $u$ is a continuous function of polynomial growth.
\end{lemma}

\begin{proof} Without loss of generality, assume that $|\phi(z)|\leq L(1+|z|)$.
Setting $w(y,z)=L(1+|y|+|z|)$,  consider the following $G$-BSDEs on $[t,T]$
\begin{align*}
\begin{split}
&U_s=\Phi(X^{t,x}_T) +\int_s^T [w( U_r, V_r)+f_0(r,X^{t,x}_r)] dr+\int_s^T [w( U_r, V_r)+g_0^{i}(r,X^{t,x}_r)] d\langle B^i,B^i\rangle_r\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ -\int_s^T  V_r dB_r- (R_T-R_s),\\
&U'_s=\Phi(X^{t,x}_T) +\int_s^T [-w( U'_r, V'_r)+f_0(r,X^{t,x}_r)] dr+\int_s^T [-w( U'_r, V'_r)+g_0^{i}(r,X^{t,x}_r)] d\langle B^i,B^i\rangle_r\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ -\int_s^T  V'_r dB_r- (R'_T-R'_s).
\end{split}
\end{align*}
By (i) of lemma \ref{lemma1} and the comparison theorem \ref{my2}, we have for each fixed $(t,x)$
\[
U'_s\leq \underline{Y}^{n,t,x}_s\leq Y^{t,x}_s\leq\bar{Y}^{n,t,x}_s \leq U_s, \ \forall s\in[t,T], \]
Recalling Proposition 4.2 in \cite{HJPS2014C}, we can find some constant $\bar{C}$ depending on $L, G,m$ and $T$ so that
\[
|U'_t|+|U_t|\leq \bar{C}(1+|x|^{q+1}),
\]
which indicates that $u$ is of   polynomial growth.

Applying (iv) of lemma \ref{lemma1} yields that for each $(t,x)\in[0,T]\times\mathbb{R}^m$,
$$\vert \underline{u}^n(t,x)-{u}(t,x)\vert\leq C_G\phi(\frac{2L}{n-L}), \ \ n>L,$$
i.e., $\underline{u}^n(t,x)$ converges to ${u}(t,x)$ uniformly in $(t,x)$. Consequently, $u$ is continuous in $(t,x)$, which ends the proof.
\end{proof}

The main result of this section is,
\begin{theorem} \label{my6}
 Let assumption \textbf{(H5)} be given. Then  $u$ is the
unique viscosity solution to the following PDE:
\begin{align}
\begin{cases}
&\partial_tu+G(H(t,x,u,D_xu,D_{x}^{2}u))+\langle
b(t,x),D_{x}u\rangle+f(t,x,u,\sigma^{\mathrm{\top}}(t,x)D_{x}u)=0,\\
&u(T,x)=\Phi(x),\ \  \ \ x\in \mathbb{R}^m,
\end{cases}
\label{feynman}
\end{align}
where
\[
\begin{array}
[c]{cl}
H_{ij}(t,x,v,p,A)=(\sigma^{\mathrm{\top}}(t,x)A\sigma(t,x))_{ij}+2\langle
p,h^{ij}(t,x)\rangle
 +2g^{i}(t,x,v,p\sigma(t,x))\mathbf{1}_{\{i=j\}}
\end{array}
\]
for any $(t,x,v,p,A)\in[0,T]\times\mathbb{R}^{m}\times\mathbb{R}
\times\mathbb{R}^{m}\times\mathbb{S}_{m}$.
\end{theorem}

For reader's convenience, we provide with the definition of  viscosity solution to equation \eqref{feynman}, see \cite{CMI}.
For every  $v\in C([0,T]\times\mathbb{R}^m)$, denote by $\mathcal{P}^{2,+}v(t,x)$ the
\textquotedblleft {parabolic superjet}\textquotedblright \ of $v$ at $(t,x)$, which refers to  the set of triples
$(a,p,X)\in \mathbb{R}\times \mathbb{R}^{m}\times \mathbb{S}_m$
such
that
\begin{align*}
v(s,y)  &  \leq v(t,x)+a(s-t)+\langle p,y-x\rangle +\frac{1}{2}\langle X(y-x),y-x\rangle+o(|s-t|+|y-x|^{2}).
\end{align*}
Similarly the
\textquotedblleft
{parabolic subjet}\textquotedblright \ of $v$ at $(t,x)$ can be defined  by $\mathcal{P}
^{2,-}v(t,x):=-\mathcal{P}^{2,+}(-v)(t,x)$.

\begin{definition}
\textup{(i)}  For $v\in C([0,T]\times \mathbb{R}^{m})$, $v$ is  a viscosity subsolution of \eqref{my6} on $[0,T]\times \mathbb{R}^{m}$, if $v(T,x)\leq \Phi(x)$ and for all $(t,x)\in(0,T)\times \mathbb{R}^{m}$,
\[
a+G({H}(t,x,{v}(t,x),p,X))+\langle b(t,x),p\rangle+{f}(t,x,{v}(t,x),p\sigma(t,x))\geq 0, \  \text{for }\ (a,p,X)\in \mathcal{P}
^{2,+}v(t,x).
\]
\textup{(ii)} A {viscosity supersolution} of equation \eqref{my6} on $[0,T]\times \mathbb{R}^{m}$ refers to function
$v\in C([0,T]\times \mathbb{R}^{m})$ with  $v(T,x)\geq \Phi(x)$ such that for each
$(t,x)\in(0,T)\times \mathbb{R}^{m}$,
\[
a+G({H}(t,x,{v}(t,x),p,X))+\langle b(t,x),p\rangle+{f}(t,x,{v}(t,x),p\sigma(t,x))\leq 0, \  \text{for }\ (a,p,X)\in \mathcal{P}
^{2,-}v(t,x).
\]

 A function $v\in C([0,T]\times \mathbb{R}^{m})$ is called a {viscosity solution} of equation \eqref{my6}
if it is simultaneously a viscosity subsolution and a viscosity supersolution of equation \eqref{my6} on $[0,T]\times \mathbb{R}^{m}$.
\end{definition}

\begin{proof}[The proof of Theorem \ref{my6}]
Since   the uniqueness of  viscosity solution to equation  \eqref{my6} is well established, c.f. \cite{PE, P10}, by the symmetry of supsolution and subsolution, we only check that $u$ is a viscosity subsolution.

Given $(t,x)\in(0, T)\times\mathbb{R}^n$ and $(a,p,X)\in\mathcal{P}^{2,+}u(t,x)$, since $\underline{u}^n$ converges to $u$ uniformly  in $(t,x)$, we get that
\begin{align*}
\lim\limits_{n\rightarrow\infty}|\underline{u}^n(t_n,x_n)-u(t,x)|=0,
\end{align*}
whenever $(t_n,x_n)\rightarrow (t,x)$.
With the help of Proposition 4.3 in \cite{CMI}, there exist  sequences
\begin{align*}
n_k\rightarrow\infty,  \ (t_k,x_k)\rightarrow(t,x),\  \text{and}\ (a_k,p_k,X_k)\rightarrow(a,p,X) ,
\end{align*}
such that
\begin{align*}
  (a_k,p_k,X_k)\in\mathcal{P}^{2,+}\underline{u}^{n_k}(t_k,x_k).
\end{align*}

From the Feynman-Kac formula in \cite{HJPS2014C},  we know $\underline{u}^n(t,x)$ is the unique viscosity solution to
\begin{align*}
\begin{cases}
&\partial_t\underline{V}+G(\underline{H}^n(t,x,\underline{V},D_x\underline{V},D_{x}^{2}\underline{V}))+\langle
b(t,x),D_{x}\underline{V}\rangle+\underline{f}_n(t,x,\underline{V},\sigma^{\top}(t,x)D_{x}\underline{V})+f_0(t,x)=0,\\
&\underline{V}(T,x)=\Phi(x),\ \  \ \ x\in \mathbb{R}^m,
\end{cases}
\end{align*}
with
\[
\begin{array}
[c]{cl}
\underline{H}^n_{ij}(t,x,v,p,A)=(\sigma^{\top}(t,x)A\sigma(t,x))_{ij}+2\langle
p,h^{ij}(t,x)\rangle
 +2\underline{g}_n^{i}(t,x,v,p\sigma(t,x))\mathbf{1}_{\{i=j\}}+2g^{i}_0(t,x)\mathbf{1}_{\{i=j\}}.
\end{array}
\]
Thus by the definition of viscosity solution we derive that
\[
a_k+G(\underline{H}^n(t_k,x_k,\underline{u}^{n_k}(t_k,x_k),p_k,X_k))+\langle b(t_k,x_k),p_k\rangle+\underline{f}_n(t_k,x_k,\underline{u}^{n_k}(t_k,x_k),p_k\sigma(t_k,x_k))+f_0(t_k,x_k)\geq 0.
\]
Recalling (i) of lemma \ref{lemma1}, we obtain that for $\varphi=f, g^{i}$
\[
\varphi(t,x,y,z)\geq \underline{\varphi}_n(t,x,y,z)+\varphi_0(t,x), \ \text{ for  all }   (t,x,y,z)\in[0,T]\times\mathbb{R}^m\times\mathbb{R}\times\mathbb{R}^{d},
\]
which implies that $H(t,x,v,p,A)\geq \underline{H}^n(t,x,v,p,A)$.
Then we derive that
\[
a_k+G({H}(t_k,x_k,\underline{u}^{n_k}(t_k,x_k),p_k,X_k))+\langle b(t_k,x_k),p_k\rangle+{f}(t_k,x_k,\underline{u}^{n_k}(t_k,x_k),p_k\sigma(t_k,x_k))\geq 0.
\]
Sending $k\rightarrow\infty$, we get
\[
a+G({H}(t,x,{u}(t,x),p,X))+\langle b(t,x),p\rangle+{f}(t,x,{u}(t,x),p\sigma(t,x))\geq 0,
\]
which is the desired result.
\end{proof}

config: arxiv


TITLE: how $\int_a ^ b |f'(x)|$ gives the length of the arc of the contour $f$ : $(f(x) : x \in [a , b])$
QUESTION [1 upvotes]: I got to know that $\int_a ^ b |f'(x)|$ gives the length of any contour. Where $f(x)$ is a piece-wise differentiable function from $[a,b]$ to $\mathbb R^2$. I was reading complex integral . Can anyone please enlighten me how it is ? I understand for the functions whose Domain and Range are $\mathbb R$.

REPLY [2 votes]: What is the length of a curve?
Let $a = x_0 < x_1 < x_2 < ... < x_n = b$ a partition of $[a,b]$
An approximation of the length can be
$Len(f) \approx \sum_{k=0}^{n-1}|f(x_i+1)-f(x_i)|$
You can see that when at the limit where the partition is delicate enough you get
$Len(f) = \sum_{k=0}^{\infty}|f(x_i+1)-f(x_i)|=\sum_{k=0}^{\infty}\Delta x_i|\frac{f(x_i+1)-f(x_i)}{\Delta x_i}|\rightarrow_{n\rightarrow \infty} = \int_a^b|f'(x)|dx$
Thats the general idea