TITLE: How to find multiplicative orders of all elements in $Z_{p^{2}}$
QUESTION [0 upvotes]: I am working on some finite fields over $Z_{p^{2}}$ and I want to compute multiplicative orders of all elements in this field. Off the top of my head, I'd say I should multiply each $p^{2}-1$ element with itself and compute modulus ${p^{2}}$ continuously till I get "1" as result. The order of that element is the number of multiplications. Is that a correct approach or not?
EDIT 1 - Suppose I want to compute multiplicative orders of all elements in $Z_{3^{2}}$:
$1 \rightarrow 1 $
$2 \rightarrow 1, 2, 4, 8, 7, 5, 1$
$3 \rightarrow 1, 3, 0$
$4 \rightarrow 1, 4, 7, 1$
$5 \rightarrow 1, 5, 7, 8, 4, 2, 1$
$6 \rightarrow 1, 6, 0$
$7 \rightarrow 1, 7, 4, 1$
$8 \rightarrow 1, 8, 1$
So, none of the elements in this field has the order equal to $p^2-1=8$. Does it mean that this field doesn't have any primitives or I am missing something?

REPLY [1 votes]: Using your notation, the question you're asking makes sense if you are referring to the multiplicative group of invertible elements, namely $G=(\mathbb{Z}_{p^2})^*$. Then it is know since Gauss that a group of this kind is cyclic and it has order $|G|=\varphi(p^2)=p(p-1)$. In your example, the order is 6 and, in fact, you have generators (2 and 5).


text:
\begin{document} \title{Flexibility Management in  Economic Dispatch with Dynamic  Automatic Generation Control}

\markboth{Submitted, {June}~2020}
{Shell \MakeLowercase{\textit{et al.}}: Bare Demo of IEEEtran.cls
for Journals}

\author{Lei Fan,~\IEEEmembership{Senior Member,~IEEE;}
		Chaoyue Zhao,~\IEEEmembership{ Member,~IEEE;}
        Guangyuan Zhang,~\IEEEmembership{Member,~IEEE;}
         Qiuhua  Huang,~\IEEEmembership{Member,~IEEE}

\thanks{Lei Fan is with the Department of Engineering Technology at University of Houston, E-mail: lfan8@uh.edu.  Chaoyue Zhao is with the Department of Industrial and Systems Engineering at University of Washington at Seattle, E-mail:cyzhao@uw.edu. Guangyuan Zhang is a senior commercial Analyst at  RWE Renewables, E-mail: zhangguangyuan1985@gmail.com. Qiuhua Huang is senior power system research engineer at Pacific Northwest National Laboratory, E-mail: qiuhua.huang@pnnl.gov.}
}

\maketitle
\vspace{-2.0cm}
\begin{abstract}
As the \lred{installation of electronically interconnected renewable energy resources grows rapidly in  power systems, system frequency maintenance \zb{and \lgreen{control} become} challenging problems \zb{to maintain the system reliability} in bulk power systems. \zb{As two of the most important frequency control actions in the control centers of independent system operators (ISOs) and utilities,} the interaction between Economic Dispatch (ED) and Automatic Generation Control (AGC) attracts more and more attention}. In this paper, we propose a robust optimization based framework to measure the system flexibility by considering the interaction between two hierarchical processes (i.e.,ED and AGC). We propose a cutting plane algorithm with the reformulation technique to obtain \lred{seven different indices} of the system. In addition, we study the impacts of several system \lred{factors (i.e., the budget of operational cost, ramping capability, and transmission line capacity)} and show numerically how these factors can \lred{influence} the system flexibility.
\end{abstract}

\begin{IEEEkeywords}
 Economic Dispatch, Automatic Generation Control (AGC), Flexibility Management, Robust Optimization, Cutting Plane Method
\end{IEEEkeywords}
\section{Nomenclature}
\addcontentsline{toc}{section}{Nomenclature}
\begin{IEEEdescription}[\IEEEusemathlabelsep\IEEEsetlabelwidth{$V_1,V_2$}]
\item [A. Sets]
\item[${\Bb}$] Set of buses.
\item[$\mathcal{{G}}$] Set of generators.
\item[$\mathcal{{G}}^b$] Set of generators at bus $b$.
\item [$\mathcal{L}$] Set of transmission lines.
\vspace{1mm}

\item[B. Parameters]

\item[$CP_n$] Penalty cost for generator $n$ in dynamic AGC constraints.
\item[$\tau$] Budget for \zb{the} total operational cost.
\item[$F_{l}$] Transmission capacity {of} transmission line $l$ {(MW)}.
\item[${\mbox{P}_n^{\mbox{\tiny max}}}$] Maximum generation amount {(MW)} of generator $n$.
\item[${\mbox{P}_{n}^{\mbox{\tiny min}}}$] Minimum generation amount {(MW)} of generator $n$.
\item[$\mbox{REG}^{\mbox{\tiny U}}_{n}$] Maximum regulation up amount of generator $n$.
\item[$\mbox{REG}^{\mbox{\tiny MinU}}$] Minimum regulation up requirement of the system.
\item[$\mbox{REG}^{\mbox{\tiny D}}_{n}$] Maximum regulation down amount of generator $n$.
\item[$\mbox{REG}^{\mbox{\tiny MinD}}$] Minimum regulation down requirement of the system.
\item[$\mbox{SR}_{n}^{\mbox{\tiny max}}$] Maximum spinning reserve amount of generator $n$.
\item[$\mbox{SR}^{\mbox{\tiny min}}$] Minimum spinning reserve requirement of the system.
\item[$\mbox{RUR}_{n}$] Maximum ramping up rate  of generator $n$.
\item[$\mbox{RDR}_{n}$] Maximum ramping down rate  of generator $n$.
\item[$\mbox{SF}_{b,l}$] Shift factors of transmission line $l$ and bus $b$.
\item[$\bar{d}_{b}$] Nominal load amount at bus $b$.
\item[$\Delta\bar{d}_{t}$] Nominal system load change amount at sub time interval $t$.
\item[$\hat{d}_{b}$] The maximum deviation amount from nominal load amount at bus $b$.
\item[$\Delta\hat{d}_{t}$] The maximum load disturbance amount from the nominal value at sub time interval $t$.
\item[$\Delta{\omega}_{t}^{\mbox{\tiny min}}$] Minimum system frequency change at sub time interval $t$.
\item[$\Delta{\omega}_{t}^{\mbox{\tiny max}}$] Maximum system frequency change at sub time interval $t$.

\vspace{1mm}

\item[C. Random Parameters]
\item[$d_{b}$] \zb{Random} load on bus $b$ (MW).
\item[$\Delta{d}_{t}$] \zb{Random} system load disturbance (MW) at sub time interval $t$.

\vspace{1mm}

\item[D. Decision Variables]

\item[${{oc}_{n}}$] {Generation} {cost} function of generator $n$.
\item[$\lambda^{up}_{b}$] The scale of upper deviation of load at bus $b$.
\item[$\lambda^{dn}_{b}$] The scale of lower deviation of load at bus $b$.
\item[$\lambda^{up}_{t}$] The scale of upper deviation of system load within one sub time interval t.
\item[$\lambda^{dn}_{t}$] The scale of lower deviation of system load within one sub time interval t.
\item[$p_{n}$] Generation amount of generator $n$ (MW).
\item[$\Delta{f}_{n,t}^{\mbox{\tiny GV}+}$] Slack variable for dynamic AGC constraints for generator $n$ at sub time interval $t$.
\item[$\Delta{f}_{n,t}^{\mbox{\tiny GV}-}$] Slack variable for dynamic AGC constraints for generator $n$ at sub time interval $t$.
\item[$reg^{\mbox{\tiny U}}_{n}$] Regulation up amount of generator $n$ (MW).
\item[$reg^{\mbox{\tiny D}}_{n}$] Regulation down amount of generator $n$ (MW).
\item[${sr_{n}}$] Spinning reserve of generator $n$ {(MW)}.
\item[$ \Delta{p}^{\mbox{\tiny GV}}_{n,t}$] The governor power generation change of generator $n$ at sub time interval $t$.
\item[$ \Delta{p}^{\mbox{\tiny M}}_{n,t}$] The prime power generation change of generator $n$ at sub time interval $t$.
\item[$\Delta{\omega}_{t}$] System frequency change.

\end{IEEEdescription}

\section{Introduction} \label{sec: Int}
As the 3D (decarbonization, digitization and decentralization) trend becomes the mainstream in the evolution of energy systems, more and more \lred{renewable energy resources} are installed in the bulk power systems.  For example, as projected by U.S. Energy Information Administration (EIA) \cite{EIA2020AEOL}, the total share of electricity generation from the renewable energy will be $38\%$ by 2050. In particular, solar energy and wind energy will contribute $17.5\%$ and $12.54\%$ of the total electricity generation by 2050 respectively. Moreover, New York state plans to reach $100\%$ \lgreen{carbon-free} by 2050 \cite{NYISO2019CP}, and California state sets \lgreen{its} $100\%$  clean electric power goal by 2045 \cite{CAISO2019CP}. The increase of these carbon-free and non-dispatchable energy resources requires the enhancement of the digital management capability of ISOs and \lred{utilities}.

In order to hedge against the variability and uncertainty of renewable energy generation and therefore achieve a high penetration of renewable energy to the power system, the concept of flexibility has been proposed and investigated, to gauge the capability of the power system in addressing the \lred{variability} of the \lred{net demand (demand net of wind and solar)} \cite{hobbs1992flx,ela2014evolution}. From the time scale's perspective, the flexibilities on planning and operations of the power system have been recently investigated. For example, the index of insufficient ramping resource expectation (IRRE) is proposed in \cite{lannoye2012evaluation} to reflect the flexibility of the power system in the generation expansion \lgreen{planning}. Operational flexibility and local flexibility for the power system operated by transmission system operators are discussed in \cite{bucher2015quantification}. From market design's perspective, ISOs developed different commodity products for the electricity market to capture the flexibility of the resources. For example, CAISO and MISO developed \lgreen{market-based} flexible ramping products that can improve the availability of  the system's ramping capacity \cite{wang2016enhancing}. ERCOT designed a fast frequency response product to maintain sufficient primary frequency control capability \lgreen{under a high} penetration of renewable energy \cite{liu2018participation}. In addition, MISO is investigating short-term reserve products to enhance the system flexibility and ensure reserve deliverability \cite{Wang2019STRP}. All these new \lred{market designs} effectively provide the pricing signal to   flexible resources in the energy market system.  On the resource level, researchers have investigated flexible generation resources such as combined-cycle power plants \cite{fan2015edge} and \cite{dai2018configuration}, pump-storage plants \cite{garcia2008stochastic}, battery energy storage \cite{ding2019rectangle} and \cite{xu2017factoring}. These complex operation models of multi-cycle or multi-stage energy resources in the electricity market not only can strengthen the capability of the grid to respond to the dynamic change of \lred{net demand}, but also can reduce the operational cost of the electricity grid \cite{guan2018unified}. In addition, the modeling approaches for aggregations of flexible resources such as virtual power plants and distributed energy resources aggregators have also been studied in \cite{babaei2019data} and \cite{chen2019aggregate}.

In this paper, we focus on the real-time flexibility of the power system, motivated by the flexibility metric framework proposed by the researchers from ISO-NE \cite{zhao2015unified}. In \cite{zhao2015unified}, the system flexibility is measured in four dimensions (i.e., time, action, uncertainty, and cost). In the time dimension, the short-term flexibility indicates the capability of the system in responding to emergencies or contingencies from minutes to hours. The long-term flexibility indicates the capability of the system in adapting the change of the generation portfolio, system topology, and regulator policy. In the action dimension, different  \lred{control schemes} (e.g., automatic generation control, economic dispatch, unit commitment, outage management, generation and transmission expansion) can be taken by system operators based on \lgreen{different response time windows} to address \lgreen{the variability} of the net demand. In the uncertainty dimension, system operators need to manage  resources to tackle \lgreen{randomnesses} such as equipment failures or forecasting errors of the net demand. Then, the cost restricts the availability of \lred{control schemes}. With these, the overall system \lred{flexibility} of the real-time economic dispatch under the uncertainty can be calculated in a systematic way.

In the current electricity market practice, ED, \lred{which usually runs every 5 minutes,}  provides generation resources with the base dispatch point and \lred{regulation reserve capacity}. However, \lred{the traditional ED process does not consider the impact of dynamic AGC, which \lgreen{is executed} every 2-6 seconds, with the objective of maintaining the system frequency. This conventional ED-AGC hierarchical model may not be able to provide sufficient system flexibility under \zb{a} high renewable penetration to track the second-to-second net demand variation  \cite{li2015miso} and \cite{Zhang2019AGCED}.} Therefore, we will study the real-time flexibility management by explicitly considering the interaction between ED and the dynamic AGC.

The contributions of our paper can be summarized as follows.
\begin{enumerate}
	\item We propose a robust optimization based framework to measure the system flexibility of the dynamic AGC constrained economic dispatch process. Then, we develop a separation framework to effectively solve the robust feasibility problem.
	\item We \lred{propose seven flexibility indices for \zb{a} systematic evaluation \zb{of the system flexibility,} and} analyze the system characteristics by using the real-time flexibility management tool to help the system operator understand the impacts of multiple system factors (e.g., budget, ramping capability, and transmission line capability) on the system flexibility.

\end{enumerate}

We organize the remaining part of this paper as follows. Section \ref{sec:MFUC} describes the  mathematical formulation of the system's real-time flexibility based on the AGC constrained economic dispatch. Section \ref{sec:CaseSt}  reports the studies of real-time flexibility in IEEE standard systems. Section \ref{sec:conclu} concludes the findings of our paper.

\section{Mathematical Formulation of  Real-Time Flexibility} \label{sec:MFUC}
{In this section, we deploy a robust optimization based framework to measure the system flexibility of the economic dispatch with the dynamic AGC. In current practice, within every 5-minute interval ED run, the system frequency is adjusted by AGC to its nominal value for every 2-6 seconds. In this paper, we consider the system dynamics within \lred{an} ED run cycle, i.e., 5 minutes. That is, the overall time horizon is set to be 5 minutes. Within the ED cycle, we consider all the AGC cycles as a set of time interval $\Tt$. Both load at each bus and system load disturbance in each time period are assumed to be random and are within \lgreen{an} undetermined variation range. That is,
\begin{align}
& {d_b} \in \mathcal{U}_b(\lambda) = [\bar{d}_{b} - \lambda^{dn}_{b}\hat{d}_b, \bar{d}_b + \lambda^{up}_{b}\hat{d}_b], \ \ \forall b \in \mathcal{B} \nonumber
\end{align}
and
\begin{align}
& \Delta{d}_{t} \in \mathcal{U}_t(\lambda) = [{\Delta\bar{d}}_{t} - \lambda^{dn}_{t}\Delta\hat{d}_t, \Delta\bar{d}_t + \lambda^{up}_{t}\Delta\hat{d}_t], \ \ \forall t \in \Tt, \nonumber
\end{align}
where $\bar{d}_{b}$ and $\hat{d}_b$ represent the nominal value and maximum deviation of the load at bus $b$, and $\Delta\bar{d}_{t}$ and $\Delta\hat{d}_t$ represent the nominal value and maximum deviation of the system load change at time period $t$. $\lambda_{\{b,t\}}^{\{up, dn\}} \in [0,1]$ represent the scales of the deviation for the corresponding variation range. Unlike the traditional robust optimization model that the size of the uncertainty set is predefined, in our model, we will obtain the largest size of each uncertainty set by deciding the scales $\lambda_{\{b,t\}}^{\{up, dn\}}$, to measure the system flexibility \lred{that} the system can accommodate.

\lred{Based on the economic dispatch with dynamic AGC model proposed by \cite{Zhang2019AGCED}, we develop a  flexibility measurement model as described in \ref{subsec:FRCTG}.}
\subsection{Formulation}\label{subsec:FRCTG}
\begin{subequations}
\begin{align}
\max \ &  \sum_{\zb{b\in \Bb}}(\hat{d}_b\lambda^{up}_{b} + \hat{d}_b\lambda^{dn}_{b}) + \sum_{t \in \Tt}(\Delta\hat{d}_t\lambda^{up}_{t} + \Delta\hat{d}_t\lambda^{dn}_{t}) \label{eqn-UC:obj}\\
s.t. &  \sum_{n \in {\mathcal{{G}}}}  oc_{n}(p_{n}) + \sum_{n \in \Gg}\sum_{t \in \Tt}\fred{\mbox{CP}_{n}\Delta{f}_{n,t}^{\mbox{\tiny GV}+}} \nonumber \\
	&  \fred{+ \sum_{\zb{n} \in \Gg}\sum_{t \in \Tt} \mbox{CP}_{n}\Delta{f}_{n,t}^{\mbox{\tiny GV}-}} \leq \tau \label{eqn-UC:budget}\\
& p_{n} + reg^{\mbox{\tiny U}}_{n} + sr_{n} \leq \mbox{P}_{n}^{\mbox{\tiny max}}, \forall n \in {\mathcal{{G}}}, \label{eqn-ED:PGUBD}\\
& \mbox{P}_{n}^{\mbox{\tiny min}} \leq p_{n} - reg^{\mbox{\tiny D}}_{n},\forall n \in {\mathcal{{G}}}, \label{eqn-ED:PGLBD}\\
& reg^{\mbox{\tiny U}}_{n} \leq \mbox{REG}^{\mbox{\tiny U}}_{n},\forall n \in {\mathcal{{G}}}, \label{eqn-ED:REGUP}\\
& reg^{\mbox{\tiny D}}_{n} \leq \mbox{REG}^{\mbox{\tiny D}}_{n},\forall n \in {\mathcal{{G}}}, \label{eqn-ED:REGDN} \\
& sr_{n} \leq \mbox{SR}_{n}^{\mbox{\tiny max}},\forall n \in \zb{\mathcal{{G}}}, \label{eqn-ED:SRUBD} \\
& \Delta{p}_{n,t+1}^{\mbox{\tiny M}} = \sum_{i \in \mathcal{G}} (\alpha_{i,n}\Delta{p}_{i,t}^{\mbox{\tiny M}} + \beta_{i,n} \Delta{p}_{i,t}^{\mbox{\tiny GV}})  + \gamma_{n} \Delta{\omega}_{t}\nonumber\\
& + \zeta_{n}\Delta{d}_{t}, \forall n \in {\mathcal{{G}}}, \forall \Delta{d}_{t} \in \mathcal{U}_t,  \forall t \in \Tt,  \label{AGC1}\\
& \Delta{\omega}_{t+1} = \sum_{i \in \mathcal{G}} (\kappa_{i}\Delta{p}_{i,t}^{\mbox{\tiny M}} + \tau_{i}\Delta{p}_{i,t}^{\mbox{\tiny GV}}) + \rho\Delta{\omega}_{t} + \eta\Delta{d}_{t}, \nonumber \\
&\forall \Delta{d}_{t} \in \mathcal{U}_t, \forall t \in \Tt, \label{AGC2}\\
& \Delta{p}_{n,t + 1}^{\mbox{\tiny GV}} - \Delta{p}_{n,t}^{\mbox{\tiny GV}} + \fred{\Delta{f}_{n,t + 1}^{\mbox{\tiny GV}+} - \Delta{f}_{n,t+1}^{\mbox{\tiny GV}-}} \nonumber \\
&= K_{n}\Delta{\omega}_{\fred{t+1}},\forall n \in {\mathcal{{G}}}, \forall t \in \Tt \label{AGC3}\\
&  \Delta{p}_{n,t + 1}^{\mbox{\tiny M}} - \Delta{p}_{n,t}^{\mbox{\tiny M}} \leq \mbox{RUR}_{n},\forall n \in {\mathcal{{G}}}, \forall t \in \Tt \label{AGC_up}\\
&  \Delta{p}_{n,t}^{\mbox{\tiny M}} - \Delta{p}_{n,t+1}^{\mbox{\tiny M}} \leq \mbox{RDR}_{n},\forall n \in {\mathcal{{G}}}, \forall t \in \Tt \label{AGC_down}\\
&  -reg^{\mbox{\tiny D}}_{n} \leq  \Delta{p}_{n,t}^{\mbox{\tiny GV}} \leq reg^{\mbox{\tiny U}}_{n}, \forall n \in {\mathcal{{G}}}, \forall t \in \Tt \label{reg}\\
&   \Delta{\omega}^{\mbox{\tiny min}}\leq \Delta{\omega}_{t} \leq \Delta{\omega}^{\mbox{\tiny max}}, \zb{\forall t \in \Tt}\label{frequency}\\
& \sum_{n\in {\mathcal{{G}}}}{p}_{n} - \sum_{b \in \Bb}{d}_{b}=0, \forall {d_b} \in \mathcal{U}_b, \label{eqn-ED-B:Pbalance} \\
&\sum_{n \in \mathcal{{G}}} sr_{n} \geq \mbox{SR}^{\mbox{\tiny min}}, \label{eqn-ED-B: spinning} \\
& \fred{\sum_{n \in \mathcal{{G}}} reg^{\mbox{\tiny U}}_{n} \geq \mbox{REG}^{\mbox{\tiny MinU}}}, \label{eqn-ED-B: Regup} \\
& \fred{\sum_{n \in \mathcal{{G}}} reg^{\mbox{\tiny D}}_{n} \geq \mbox{REG}^{\mbox{\tiny MinD}}},\label{eqn-ED-B: Regdn} \\
& - F_{l} \leq \sum_{b \in \mathcal{B}} \mbox{SF}_{b,l}(\sum_{n \in \mathcal{G}^{b}}p_{n} - {d}_{b}) \leq F_{l}, \forall {d_b} \in \mathcal{U}_b, \label{transmissioncap}\\
&  p_{n},  reg^{\mbox{\tiny U}}_{n},  reg^{\mbox{\tiny D}}_{n}, {sr_{n}}, \zb{\Delta{f}_{n,t}^{\mbox{\tiny GV}+}, \Delta{f}_{n,t}^{\mbox{\tiny GV}-}} \geq 0, \nonumber\\
& \Delta{p}^{\mbox{\tiny GV}}_{n,t},\Delta{p}^{\mbox{\tiny M}}_{n,t},\Delta{\omega}_{t} \ \mbox{free},\\
& \lambda_b^{up}, \lambda_b^{dn},  \lambda_t^{up}, \lambda_t^{dn} \in [0,1], \zb{\forall b \in \Bb, \forall n \in {\mathcal{{G}}}}, \forall t \in \Tt, \nonumber
\end{align}
\end{subequations}
where the objective  function is to maximize the variation range of the uncertainty. Constraints \eqref{eqn-UC:budget} represent the budget constraint, which indicates that the total fuel cost and the penalty cost should not exceed a budget $\tau$. Constraints (\ref{eqn-ED:PGUBD}) and (\ref{eqn-ED:PGLBD}) represent {the generation} limits of {traditional} thermal units that \lgreen{take} account of generation output, regulation and spinning reserves. Constraints (\ref{eqn-ED:REGUP}), (\ref{eqn-ED:REGDN}), and (\ref{eqn-ED:SRUBD}) represent the capacities for providing regulation up, regulation down, and spinning reserve services respectively. Constraints \eqref{AGC1} - \eqref{AGC3} represent AGC dynamic system constraints, i.,e., the transformation of state vectors $\Delta{p}_{n,t}^{\mbox{\tiny M}}$, $\Delta{p}_{n,t}^{\mbox{\tiny GV}}$, $\Delta{\omega}_t$ from $t$ to $t+1$ given any demand disturbance $\Delta{d}_t$, where the matrix components $\alpha, \beta, \gamma, \zeta, \kappa, \tau, \rho, \eta$ can be calculated by numerous methods such as zero-order hold method \cite{Zhang2019AGCED}. Constraints \eqref{AGC_up} and \eqref{AGC_down} \lred{restrict} ramping up and ramping down limits. Constraints \eqref{reg} indicate that the governor generation change should not exceed the regulation service reserved, and constraints \eqref{frequency} restrict the limit of system frequency change.
In addition, the power balance constraints are described in (\ref{eqn-ED-B:Pbalance}), which should be held for any load realization within the uncertainty set; constraints
(\ref{eqn-ED-B: spinning})-\eqref{eqn-ED-B: Regdn} describe the overall spinning reserve, regulation up, and regulation down requirements respectively, and constraints \eqref{transmissioncap} represent the transmission capacity constraints.

\subsection{Solution Methodology}
First, for notation brevity, we use matrices and vectors to represent constraints and variables, and rewrite the above model in \lgreen{an} abstract compact form (denoted as ACF):
\begin{subequations}
\begin{align}
(\mbox{ACF})\ \  & \max \ \zb{a^T} \lambda \label{primal_abs_obj}\\
\hspace{-0.1cm}s.t.\ \ &  A_1 x \leq b_{1}, \label{primal_cons1} \\
& A_2 x = H_2 d, \ \ \zb{\forall d_b \in \mathcal{U}_b(\lambda), \forall b \in \Bb}, \label{primal_cons6}\\
&A_3 x \leq H_3 d, \ \ \zb{\forall d_b \in \mathcal{U}_b(\lambda), \forall b \in \Bb}, \label{primal_cons7}\\
&\zb{A_{4}y = H_{4} \Delta d, \ \ \Delta d_t \in \mathcal{U}_t(\lambda), \zb{\forall t \in \Tt},} \label{primal_cons2}\\
&\zb{A_{5}y = b_{5},  }\label{primal_cons3}\\
&\zb{A_{6} y \leq b_{6}, }\label{primal_cons4}\\
&\zb{A_7x + A_{8} y \leq b_{7},} \label{primal_cons5}
\end{align}
\end{subequations}
where \zb{$x = (p, reg^{\mbox{\tiny U}}, reg^{\mbox{\tiny D}}, sr)$, $y = (\Delta p^{\mbox{\tiny M}}$, $\Delta{p}^{\mbox{\tiny GV}}$, $\Delta{\omega}$, $\Delta{f}^{\mbox{\tiny GV}+}$, $\Delta{f}^{\mbox{\tiny GV}-})$, $d = (d_1, \cdots,$ $d_b, \cdots)_{b \in \Bb}$, and $\Delta d = (\Delta d_1$, $\cdots$, $\Delta d_t$, $\cdots)_{t \in \Tt}$}; objective \eqref{primal_abs_obj} \lgreen{represents} \eqref{eqn-UC:obj}; constraint \eqref{primal_cons1} represents \eqref{eqn-ED:PGUBD} - \eqref{eqn-ED:SRUBD} and \eqref{eqn-ED-B: spinning} - \eqref{eqn-ED-B: Regdn}; constraint \eqref{primal_cons6} represents \eqref{eqn-ED-B:Pbalance}; constraint \eqref{primal_cons7} represents \eqref{transmissioncap}; constraint \eqref{primal_cons2} represents \eqref{AGC1} and \eqref{AGC2}; constraint \eqref{primal_cons3} represents \eqref{AGC3}; constraint \eqref{primal_cons4} represents \eqref{AGC_up}, \eqref{AGC_down}, and \eqref{frequency}; constraint \eqref{primal_cons5} represents \eqref{eqn-UC:budget} and \eqref{reg}.

We deploy Benders' decomposition framework to solve the problem. Since the problem is to find the largest deviation of the uncertainty set that the system can accommodate, i.e., the largest value of $\lambda$ without making the constraints \eqref{primal_cons1}-\eqref{primal_cons5} infeasible, therefore, only feasibility cuts are needed.
\subsubsection{Master Problem and Subproblem}
We first decompose problem ACF into a master problem   \lred{(denoted as MAP)} and a subproblem.
\begin{subequations}
\begin{align}
(\mbox{MAP})~~ \max \ & \zb{a^T} \lambda \label{primal_MAP_obj}\\
s.t.~~&  g(\lambda) \leq 0 \label{primal_MAP_cons1} \\
& \textbf{0} \leq \lambda \leq \zb{\textbf{1}}
\end{align}
\end{subequations}
Here, the feasibility cuts are represented in \eqref{primal_MAP_cons1}, \zb{and \textbf{0}, $\textbf{1}$ represent vectors with all components 0 and 1 respectively}.

To generate the feasibility cuts, we first describe the feasibility check problem \lred{ (denoted as FEA)} as follows:
\begin{subequations}
\begin{align}
(\mbox{FEA}) & \max_{\zb{d_b} \in \mathcal{U}_b, \Delta \zb{d_t} \in \mathcal{U}_t}  \min_{x,y,s}  \ \zb{\textbf{1}}^T s \label{Fea_abs_obj}\\
s.t. \ \  & \fred{A_1} x - s_1\leq b_{1}, \label{Fea_cons1} \\
  &A_2 x +s_{2}^+ - s_{2}^- = H_2 d, \label{Fea_cons6}\\
 &A_3 x  - \zb{s_3} \leq H_3 d, \label{Fea_cons7}\\
 &\zb{A_{4}}y + s_{4}^+ - s_{4}^-  = H_{4} \Delta d, \label{Fea_cons2}\\
 &\fred{\zb{A_{5}}y + s_{5}^+ - s_{5}^- = b_{5}}, \label{Fea_cons3}\\
 &\zb{A_{6}} \zb{y} - s_{6} \leq b_{6}, \label{Fea_cons4}\\
 &A_7x + {\zb{A_{8}}} \zb{y} - s_{7} \leq b_{7}. \label{Fea_cons5}
 \end{align}
 \end{subequations}
If $\lambda$ is feasible, then the optimal value of (FEA) will be $0$.

\subsubsection{Reformulation of Subproblem}
 Now we take the dual of the inner minimization of subproblem (FEA) and combine the dual problem with the outer maximization problem, then we can get the following formulation:
\begin{subequations}
\begin{align}
(\mbox{DFEA}) &\max_{\zb{d_b} \in \mathcal{U}_b, \Delta \zb{d_t} \in \mathcal{U}_t, \mu}  \   \zb{b_1^T \mu_1  + d^T H_2^T \mu_2  + d^T H_3^T \mu_3}   \nonumber \\
& \zb{+  \Delta d^T H_{4}^T \mu_{4} + b_{5}^T \mu_{5} + b_{6}^T \mu_{6} + b_{7}^T \mu_{7}}\\
s.t.\ &  \zb{A_1^T \mu_1 + A_2^T \mu_2 + A_3^T \mu_3 + A_{7}^T \mu_{7}}  \leq 0, \\
& \zb{A_{4}^T}\mu_{4} + \zb{A_{5}^T}\mu_{5} + \zb{A_{6}^T}\mu_{6} + \zb{A_{8}^T}\mu_{7} = 0, \\
& \fred{ -\zb{\textbf{1}} \leq \mu_1, \mu_3, \mu_6, \mu_7 \leq \textbf{0},  } \\
& \fred{ -\zb{\textbf{1}} \leq \mu_2, \mu_{4}, \mu_{5} \leq \zb{\textbf{1}}},
\end{align}
\end{subequations}
where $\mu_1$, $\mu_2$, $\mu_3$, $\mu_{4}$, $\mu_{5}$, $\mu_{6}$, $\mu_{7}$ are dual variables for constraints \eqref{Fea_cons1} - \eqref{Fea_cons5} respectively.

\fred{In the above formulation (DFEA), we have  bilinear terms \zb{$d^T H_2^T \mu_2$, $d^T H_3^T \mu_3$, and $\Delta d^T H_{4}^T \mu_{4}$}. We will deal with \zb{$d^T H_3^T \mu_3$} first. \zb{Let $N_i$ represent the dimension of $\mu_i$ and $\lambda^{*}$ represent} the optimal solution of problem (MAP). Based on the property of $d$, \zb{$\forall b \in \Bb$,} we can rewrite it as \zb{$d_b = \bar{d}_b + z^{+}_b \lambda^{*,up}_b \hat{d}_b -  z^{-}_b \lambda^{*,dn}_b \hat{d}_b$}. Here we introduce two binary variables $\zb{z^{+}_b}$ and $z^{-}_b$ to indicate the deviation direction. Note that the variables  $\lambda^{*,up}$ and $\lambda^{*,dn}$ have been fixed for DFEA problem.  Therefore, we can replace the bilinear term $d^T H_3^T \mu_3$ as follows:}
\begin{subequations}
\begin{align}
& d^T H_3^T \mu_3 = \sum_{b\in \Bb}\sum_{i=1}^{N_3} d_b \lgreen{H_{3,i,b}}\mu_{3,i}\nonumber\\
& = \sum_{b \in \Bb}\sum_{i=1}^{N_3}(\bar{d}_b \lgreen{H_{3,i,b}} \mu_{3,i} + \lambda^{*,up}_b\hat{d}_b \lgreen{H_{3,i,b}} \zb{z^{+}_b}\mu_{3,i} \nonumber\\
&- \lambda^{*,dn}_b\hat{d}_b \lgreen{H_{3,i,b}} \zb{z^{-}_b}{\mu}_{3,i}) \label{Refm_dmu}\\
& \zb{z^{+}_b} + \zb{z^{-}_b} = 1, z_b^{+}, z_b^{-} \in \{0,1\} \label{Refm_unique}
\end{align}
\end{subequations}
\fred{
In (\ref{Refm_dmu}), we have two bilinear items $\zb{z^{+}_b}\mu_{3,i}$ and $\zb{z^{-}_b}\mu_{3,i}$, which can be further linearized by introducing auxiliary variables $\mu_{3,b,i}^{+}$ and $\mu_{3,b,i}^{-}$.
By following the approach indicated in \cite{Fan2014MinMaxRegret}, we have the following reformulation which is equivalent to (\ref{Refm_dmu}):}
\zb{
\begin{subequations}
\begin{align}
&  d^T H_3^T \mu_3 = \sum_{b \in \Bb}\sum_{i=1}^{N_3}(\bar{d}_b \lgreen{H_{3,i,b}}\mu_{3,i} + \nonumber\\
&\lambda^{*,up}_b\hat{d}_b \lgreen{H_{3,i,b}} \mu_{3,b,i}^{+} - \lambda^{*,dn}_b\hat{d}_b \lgreen{H_{3,i,b}} {\mu}_{3,b,i}^{-}) \label{Refm_dmuR}\\
&-\zb{z^{+}_b} \leq \mu_{3,b,i}^{+}, \ \mu_{3,i} \leq \mu_{3,b,i}^{+} \leq 1 - \zb{z^{+}_b} + \mu_{3,i},  \\
&-\zb{z^{-}_b} \leq \mu_{3,b,i}^{-},\ \mu_{3,i} \leq \mu_{3,b,i}^{-} \leq 1 - \zb{z^{-}_b} + \mu_{3,i}, \\
& z_b^+ + z_b^- = 1, z_b^{+}, z_b^{-} \in \{0,1\},  \\
&-1 \leq \mu_{3,b,i}^{+}, \mu_{3,b,i}^{-} \leq 0, \forall b \in \Bb, \forall i = 1, \cdots, N_3. \label{last}
\end{align}
\end{subequations}
}
For \zb{$d^T H_2^T \mu_2$ and $\Delta d^T H_{4}^T \mu_{4}$}, we can use a similar approach. \zb{We will use $d^T H_2^T \mu_2$ as an example and the other follows the same approach. Since $\mu_2 \in [-\textbf{1},\textbf{1}]$, we can replace it with $\mu_{2}^n - \mu_{2}^p$ and $- \textbf{1} \leq \mu_{2}^n \leq \textbf{0}$ and $-\textbf{1}\leq \mu_{2}^p \leq \lgreen{\textbf{0}} $. Then $d^T H_2^T\mu_2 = d^T H_2^T\mu_{2}^n - d^T H_2^T\mu_{2}^p$. For $d^T H_2^T\mu_{2}^n$ and $d^T H_2^T\mu_{2}^p$}, we will follow the same procedure of \eqref{Refm_dmuR}-\eqref{last} to linearize them. Therefore, we can reformulate the (DEFA) problem as follows:
\zb{
\begin{subequations}
\begin{align}
(\mbox{RDFEA}) \ &\max_{d_b, \Delta d_t, \mu}  \   b_1^T \mu_1   +\sum_{b \in \Bb}\sum_{i=1}^{N_2} \{ \bar{d}_b \lgreen{ H_{2,i,b}} (\mu_{2,i}^n-\mu_{2,i}^p)+ \nonumber\\
& \hspace{-1.6cm}\lambda^{*,up}_b\hat{d}_b\lgreen{ H_{2,i,b}} (\mu_{2,b,i}^{n,+}-\mu_{2,b,i}^{p,+}) - \lambda^{*,dn}_b\hat{d}_b \lgreen{ H_{2,i,b}} ({\mu}_{2,b,i}^{n,-}-{\mu}_{2,b,i}^{p,-})\} \nonumber\\
&\hspace{-1.6cm}  + \sum_{b \in \Bb}\sum_{i=1}^{N_3}(\bar{d}_b \lgreen{ H_{3,i,b}} \mu_{3,i} + \lambda^{*,up}_b\hat{d}_b \lgreen{ H_{3,i,b}}\mu_{3,b,i}^{+} \nonumber\\
&\hspace{-1.6cm} - \lambda^{*,dn}_b\hat{d}_b \lgreen{ H_{3,i,b}} {\mu}_{3,b,i}^{-}) + \sum_{t \in \Tt}\sum_{i=1}^{N_4} \{ \Delta \bar{d}_t \lgreen{ H_{4,i,t}}(\mu_{4,i}^n-\mu_{4,i}^p)  \nonumber\\
& \hspace{-1.6cm}+\lambda^{*,up}_t \Delta\hat{d}_t \lgreen{ H_{4,i,t}}(\mu_{4,t,i}^{n,+}- \mu_{4,t,i}^{p,+}) \nonumber\\
& \hspace{-1.6cm}- \lambda^{*,dn}_b \lgreen{\Delta\hat{d}_t} \lgreen{ H_{4,i,t}} ({\mu}_{4,t,i}^{n,-}-{\mu}_{4,t,i}^{p,-})\} + b_{5} \mu_{5}  + b_{6} \mu_{6} + b_{7} \mu_{7} \nonumber\\
& \hspace{-1cm}  s.t. \ \ \ \lgreen {A_1^{T} \mu_1 + A_2^{T} \mu_2 + A_3^{T} \mu_3 + A_{7}^{T} \mu_{7}  \leq 0}, \\
& \hspace{-0.8 cm} \lgreen{ \zb{A_{4}^{T}}\mu_{4} + \zb{A_{5}^{T}}\mu_{5} + \zb{A_{6}^{T}}\mu_{6} + \zb{A_{8}^{T}}\mu_{7} = 0},\\
&\hspace{-0.8cm}  -\zb{z^{+}_b} \leq \mu_{2,b,i}^{n,+}, \ \mu_{2,i}^n \leq \mu_{2,b,i}^{n,+} \leq 1 - \zb{z^{+}_b} + \mu_{2,i}^n,\\
&\hspace{-0.8cm}  -\zb{z^{-}_b} \leq \mu_{2,b,i}^{n,-},\  \mu_{2,i}^n \leq \mu_{2,b,i}^{n,-} \leq 1 - \zb{z^{-}_b} + \mu_{2,i}^n,\\
&\hspace{-0.8cm}  -\zb{z^{+}_b} \leq \mu_{2,b,i}^{p,+},\  \mu_{2,i}^p \leq \mu_{2,b,i}^{p,+} \leq 1 - \zb{z^{+}_b} + \mu_{2,i}^p,\\
&\hspace{-0.8cm}  -\zb{z^{-}_b} \leq \mu_{2,b,i}^{p,-}, \ \mu_{2,i}^p \leq \mu_{2,b,i}^{p,-} \leq 1 - \zb{z^{-}_b} + \mu_{2,i}^p,\\
&\hspace{-0.8cm} \forall b \in \Bb, \forall i = 1, \cdots, N_2, \nonumber\\
&\hspace{-0.8cm}  -\zb{z^{+}_b} \leq \mu_{3,b,i}^{+},\  \mu_{3,i} \leq \mu_{3,b,i}^{+} \leq 1 - \zb{z^{+}_b} + \mu_{3,i},\\
&\hspace{-0.8cm}  -\zb{z^{-}_b} \leq \mu_{3,b,i}^{-},\  \mu_{3,i} \leq \mu_{3,b,i}^{-} \leq 1 - \zb{z^{-}_b} + \mu_{3,i},\\
&\hspace{-0.8cm} \forall b \in \Bb, \forall i = 1, \cdots, N_3, \nonumber\\
&\hspace{-0.8cm}  -\zb{z^{+}_t} \leq \mu_{4,t,i}^{n,+},\  \mu_{4,i}^n \leq \mu_{4,t,i}^{n,+} \leq 1 - \zb{z^{+}_t} + \mu_{4,i}^n,\\
&\hspace{-0.8cm}  -\zb{z^{-}_t} \leq \mu_{4,t,i}^{n,-},\  \mu_{4,i}^n \leq \mu_{4,t,i}^{n,-} \leq 1 - \zb{z^{-}_t} + \mu_{4,i}^n,\\
&\hspace{-0.8cm}  -\zb{z^{+}_t} \leq \mu_{4,t,i}^{p,+},\  \mu_{4,i}^p \leq \mu_{4,t,i}^{p,+} \leq 1 - \zb{z^{+}_t} + \mu_{4,i}^p,\\
&\hspace{-0.8cm}  -\zb{z^{-}_t} \leq \mu_{4,t,i}^{p,-},\ \mu_{4,i}^p \leq \mu_{4,t,i}^{p,-} \leq 1 - \zb{z^{-}_t} + \mu_{4,i}^p,\\
&\hspace{-0.8cm} \forall t \in \Tt, \forall i = 1, \cdots, N_4, \nonumber\\
&\hspace{-0.8cm}   z_b^+ + z_b^- = 1, \ \ z_t^+ + z_t^- = 1, \forall b \in \Bb, \forall t \in \Tt, \\
&\hspace{-0.8cm}  -\textbf{1} \leq\mu_{2}^{n,+}, \mu_{2}^{n,-}, \mu_{2}^{p,+}, \mu_{2}^{p,-}  \leq \textbf{0}, \\
&\hspace{-0.8cm} -\textbf{1}\leq  \mu_{3}^{+}, \mu_{3}^{-} \leq \textbf{0}, \nonumber\\
&\hspace{-0.8cm}  -\textbf{1} \leq\mu_{4}^{n,+}, \mu_{4}^{n,-} , \mu_{4}^{p,+},\mu_{4}^{p,-} \leq \textbf{0}, \\
&\hspace{-0.8cm}   -\textbf{1} \leq \mu_1, \mu_3, \mu_6, \mu_7 \leq \textbf{0},  \\
&\hspace{-0.8cm}   -\textbf{1} \leq \mu_2, \mu_{4}, \mu_{5} \leq \textbf{1}.
\end{align}
\end{subequations}
}
If the objective value $\eta^{*}$ of subproblem (RDFEA) is greater than $0$, we can obtain the following feasibility cut:

\zb{
\begin{align}
& \eta(\lambda) = \nonumber \\
&\  \sum_{b \in \Bb}(\sum_{i=1}^{N_2}\hat{d}_b \lgreen{H_{2,i,b}}(\mu_{2,b,i}^{n+,*}-\mu_{2,b,i}^{p+,*}) +  \sum_{i=1}^{N_3}\hat{d}_b\lgreen{H_{3,i,b}}\mu_{3,b,i}^{+,*})\lambda^{up}_b   \nonumber\\
& -  \sum_{b \in \Bb}(\sum_{i=1}^{N_2}\hat{d}_b \lgreen{H_{2,i,b}}(\mu_{2,b,i}^{n-,*}-\mu_{2,b,i}^{p-,*}) +  \sum_{i=1}^{N_3}\hat{d}_b \lgreen{H_{3,i,b}}\mu_{3,b,i}^{-,*})\lambda^{dn}_b   \nonumber\\
&+\sum_{t \in \Tt}\sum_{i=1}^{N_4} \Delta\hat{d}_t \lgreen{H_{4,i,t}}(\mu_{4,t,i}^{n+,*}-\mu_{4,t,i}^{p+,*})\lambda^{up}_t\nonumber\\
& - \sum_{t \in \Tt}\sum_{i=1}^{N_4} \lgreen{\Delta\hat{d}_t}\lgreen{H_{4,i,t}} ({\mu}_{4,t,i}^{n-,*}-{\mu}_{4,t,i}^{p-,*})\lambda^{dn}_t \nonumber\\
& +b_1^T \mu_1^{*}   +\sum_{b \in \Bb}\sum_{i=1}^{N_2} \bar{d}_b \lgreen{H_{2,i,b}}(\mu_{2,i}^{n*}-\mu_{2,i}^{p*})\nonumber\\
& + \sum_{b \in \Bb}\sum_{i=1}^{N_3}\bar{d}_b \lgreen{H_{3,i,b}}\mu_{3,i}^{*} + \sum_{t \in \Tt}\sum_{i=1}^{N_4}  \Delta\bar{d}_t\lgreen{H_{4,i,t}}(\mu_{4,i}^{n*}-\mu_{4,i}^{p*})\nonumber\\
& + b_{5} \mu_{5}^{*}  + b_{6} \mu_{6}^{*} + b_{7} \mu_{7}^{*} \leq 0 \label{feaCUT}
\end{align}}
\subsubsection{Algorithm Framework}
We summarize our Benders' decomposition based framework as follows:
\begin{enumerate}
	\item Solve the master problem (MAP) and obtain the optimal solution $\lambda^*$.
	\item Test the feasibility of the subproblem by solving dual reformulated DFEA problem with $\lambda^*$.
	\item If the optimal value $\eta^{*}$ of DFEA $ > 0$, then update the master problem by adding feasibility cut \eqref{feaCUT} and go to Step 2. Otherwise, terminate and return the optimal solution and objective value of master problem.
\end{enumerate}

\begin{figure}[tph]
	\begin{center}

		\usetikzlibrary{shapes,arrows}

		\tikzstyle{decision} = [diamond, draw,,outer color=blue!60,inner color=white,very thick, fill=white!20, text
		width=6em, text badly centered, node distance=4cm, inner sep=0pt]
		\tikzstyle{smalldecision} = [diamond, draw, fill=white!20, text
		width=4em, text badly centered, node distance=4cm, inner sep=0pt]
		\tikzstyle{block} = [rectangle, draw,outer color=brown!60,inner color=white,very thick, fill=white!20, rounded
		corners, minimum size=1em] \tikzstyle{bigblock} = [rectangle, draw,
		fill=white!20, rounded corners, text width=14em, minimum height=3em]
		\tikzstyle{midblock} = [rectangle, draw,,outer color= green!60,inner color=white,very thick, fill=white!20, rounded
		corners, text width=6em, minimum height=4em] \tikzstyle{line} =
		[draw, ->, thick, color=black] \tikzstyle{cloud} = [draw,
		ellipse,fill=black, minimum size=1em]
		\tikzstyle{point}=[circle,inner sep=0pt,minimum size=2pt,fill=black]
		\tikzstyle{every node}=[font=\large]
		\tikzset{
			every node/.style={scale=0.8}
		}

		\begin{tikzpicture}[scale=1.2, node distance = 2cm, inner sep=6pt, text badly centered,auto]

		\node [block] (init) at (2.5,10.5) {{Initialization}};
		\node [midblock] (Smaster) at (2.5,9) {Solve $\mbox{MAP}$ and {obtain $\lambda^*$}};
		\node [midblock] (Gsub) at (2.5,7.75) {Formulate $\mbox{RDFEA}$ and {obtain $\eta^{*}$}};
		\node [decision] (Cut) at (2.5, 6) {$\eta^{*} > 0$ ?};
		\node [midblock] (Update) at (0.2, 7.75){Add $\eta <= 0$ into $\mbox{MAP}$};
		\node [block] (stop) at (2.5, 4.5){Terminate and output};

		\path [line] (init) -- (Smaster);
		\path [line] (Smaster) -- (Gsub);
		\path [line] (Gsub) -- (Cut);
		\path [line] (Update) |- (Smaster);
		\path [line] (Cut) -| node [near start] {Yes} (Update);
		\path [line] (Cut) -- node {No}(stop);

		\end{tikzpicture}
	\end{center}
	\caption{{{Flowchart} of Decomposition Algorithm}}\label{ALG:BendersChart}
\end{figure}

\section{Case Study}\label{sec:CaseSt}
In this section, we test the performance of the proposed real-time flexibility management model with the IEEE 118-bus system (online \cite{IIT_Systems}). These cases are tested by using Julia \cite{lubin2015computing} and CPLEX 12.8 \cite{IBMCPLEX} on Intel Xeon Silver 4216 CPU and 128 G memory.

To get the insights \lgreen{on} how the system operators can improve their system flexibility by using limited resources, we investigate on which factors will have impacts on the system flexibility. More specifically, we study key factors including the budget of operational cost, ramping capability, and transmission line capacity. In order to obtain more information on the system's capability in handling uncertainty, we analyze seven different system flexibility indicators: \lred{total flexibility (TF), economic dispatch flexibility (EDF), AGC flexibility (AGCF), economic dispatch upward flexibility (EDUPF), economic dispatch downward flexibility (EDDNF), AGC upward flexibility (AGCUPF), and AGC downward flexibility (AGCDNF)} based on the simulation results. The definitions of these flexibility indicators are shown in (\ref{eqn-TF}) - (\ref{eqn-AGCDNF}).
\begin{subequations}
  \begin{align}
 & \mbox{TF}  = \sum_{b}(\hat{d}_b\lambda^{up}_{b} + \hat{d}_b\lambda^{dn}_{b}) + \sum_{t \in \Tt}(\Delta\hat{d}_t\lambda^{up}_{t} + \Delta\hat{d}_t\lambda^{dn}_{t}) \label{eqn-TF}\\
 & \mbox{EDF}  = \sum_{b}(\hat{d}_b\lambda^{up}_{b} + \hat{d}_b\lambda^{dn}_{b})  \label{eqn-EDF}\\
 & \mbox{AGCF}  =  \sum_{t \in \Tt}(\Delta\hat{d}_t\lambda^{up}_{t} + \Delta\hat{d}_t\lambda^{dn}_{t}) \label{eqn-AGCF}\\
 & \mbox{EDUPF}  = \sum_{b}(\hat{d}_b\lambda^{up}_{b})  \label{eqn-EDUPF}\\
 & \mbox{EDDNF}  =  \sum_{b}(\hat{d}_b\lambda^{dn}_{b})  \label{eqn-EDDNF}\\
 & \mbox{AGCUPF}  =  \sum_{t \in \Tt}(\Delta\hat{d}_t\lambda^{up}_{t} ) \label{eqn-AGCUPF}\\
 & \mbox{AGCDNF}  =  \sum_{t \in \Tt}(\Delta\hat{d}_t\lambda^{dn}_{t}) \label{eqn-AGCDNF}
  \end{align}
  \end{subequations}
The TF can gauge both static and dynamic system capability to handle the uncertainty caused by forecasting errors of load and renewable energy resources. The EDF reflects the static capability of the system. The AGCF indicates the capability of the system to respond to the disturbance in real time. The EDUPF and EDDNF can explain up and down static flexibility of the system.  The \lred{AGCUPF} and AGCDNF can describe the up and down dynamic flexibility of the system, respectively.

In the nominal case, there are 30 units online and the system load is  \lgreen{around} $86.3\%$ of the total capacity of online units. We assume the uncertain range of the net load at each bus is $[-15\%, 15\%]$ and the system load disturbance is within the range $[-1\%, 1\%]$.

\subsection{Impact of Operational Cost}
In this subsection, we first investigate the impact of the operational cost budget by increasing it gradually (denoted from B0 to B22 in Column 1, Table \ref{tab:SFUDB}).  We set the budget of the base case (denoted as B0) as the optimal operational cost of the economic dispatch  under the nominal foretasted net load without the dynamic AGC. Then we incrementally scale up the budget by using the scale factor (SF) as shown in Column 2 in Table \ref{tab:SFUDB}.  For example, the budget of B22 is four times of the budget of B0.  All other system parameters are the same in these 23 cases. We report the simulation results of the proposed flexibility indicators under different budgets in Table \ref{tab:SFUDB}.

\begin{table}[h]
	\huge
	\centering
	\caption{{System Flexibility under Different Budgets} }
	\label{tab:SFUDB}
	\begin{adjustbox}{max width=0.48\textwidth}
		\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
			\hline
	Budget&	SF & TF      & EDF     & AGCF  & EDUPF  & EDDNF  & AGCUPF & AGCDNF \\ \hline
	B0	&	1.00        & 115.298 & 115.298 & 0.000 & 47.626 & 67.672 & 0.000  & 0.000  \\ \hline
	B1	&	1.01        & 120.678 & 120.678 & 0.000 & 53.006 & 67.672 & 0.000  & 0.000  \\ \hline
	B2	&	1.02        & 125.532 & 125.532 & 0.000 & 57.861 & 67.672 & 0.000  & 0.000  \\ \hline
	B3	&	1.03        & 129.411 & 129.411 & 0.000 & 61.739 & 67.672 & 0.000  & 0.000  \\ \hline
	B4	&	1.04        & 129.520 & 129.476 & 0.044 & 61.804 & 67.672 & 0.026  & 0.018  \\ \hline
	B5	&	1.05        & 129.550 & 129.476 & 0.074 & 61.804 & 67.672 & 0.048  & 0.026  \\ \hline
	B6	&	1.06        & 129.579 & 129.476 & 0.103 & 61.804 & 67.672 & 0.078  & 0.026  \\ \hline
	B7	&	1.07        & 129.609 & 129.476 & 0.133 & 61.804 & 67.672 & 0.059  & 0.074  \\ \hline
	B8	&	1.08        & 129.638 & 129.476 & 0.162 & 61.804 & 67.672 & 0.086  & 0.076  \\ \hline
	B9	&	1.09        & 129.668 & 129.476 & 0.192 & 61.804 & 67.672 & 0.093  & 0.099  \\ \hline
	B10	&	1.10        & 129.697 & 129.476 & 0.221 & 61.804 & 67.672 & 0.109  & 0.112  \\ \hline
	B11	&	1.20        & 129.993 & 129.476 & 0.517 & 61.804 & 67.672 & 0.290  & 0.227  \\ \hline
	B12	&	1.30        & 130.288 & 129.476 & 0.812 & 61.804 & 67.672 & 0.400  & 0.412  \\ \hline
	B13	&	1.40       & 130.583 & 129.476 & 1.107 & 61.804 & 67.672 & 0.578  & 0.529  \\ \hline
	B14	&	1.50        & 130.878 & 129.476 & 1.402 & 61.804 & 67.672 & 0.730  & 0.672  \\ \hline
	B15	&	1.60        & 131.174 & 129.476 & 1.698 & 61.804 & 67.672 & 0.682  & 1.015  \\ \hline
	B16	&	1.70        & 131.469 & 129.476 & 1.993 & 61.804 & 67.672 & 0.844  & 1.149  \\ \hline
	B17	&	1.80        & 131.764 & 129.476 & 2.288 & 61.804 & 67.672 & 1.133  & 1.156  \\ \hline
	B18	&	1.90        & 131.976 & 129.476 & 2.500 & 61.804 & 67.672 & 1.250  & 1.250  \\ \hline
	B19	&	2.00        & 131.976 & 129.476 & 2.500 & 61.804 & 67.672 & 1.250  & 1.250  \\ \hline
	B20	&	2.50        & 131.976 & 129.476 & 2.500 & 61.804 & 67.672 & 1.250  & 1.250  \\ \hline
	B21	&	3.00        & 131.976 & 129.476 & 2.500 & 61.804 & 67.672 & 1.250  & 1.250  \\ \hline
	B22	&	4.00        & 131.976 & 129.476 & 2.500 & 61.804 & 67.672 & 1.250  & 1.250  \\ \hline
		\end{tabular}
	\end{adjustbox}
\end{table}

From Columns 1-5 in Table \ref{tab:SFUDB}, we can observe that three flexibility indicators (i.e., TF, EDF, and \lred{AGCF}) are non-decreasing as the budget increases.
On one hand, for EDF, when the budget is small and gradually increases by $1\%$ from B0 to B3, the EDF increases by around $4\%$ for each level of the budget. Once \lgreen{the} scale factor (SF) reaches $1.04$, EDF remains as a constant number. On the other hand, when the budget remains small, i.e., B0-B3, AGCF remains 0. But when the budget becomes larger, AGCF increases and becomes a constant \lgreen{when a relatively large budget (i.e., B18 and above) is available.} That is because that when the budget is small, EDF will be allocated to all operational cost budget as it has more weight than AGCF in the objective function. But when the system has more budget for operational cost,  which is sufficient for the operational cost of the economic dispatch, it will start to allocate budget to the dynamic AGC part. Therefore, with more budget available, AGCF will incur more operational \lgreen{costs}. But when the total budget is sufficient for both ED and AGC costs, both EDF and AGCF will not increase. In addition, the increment of AGCF will result in the increment of TF. We show the trends of indicators TF, EDF, and AGCF in Fig. \ref{FIG:sfB}.

\begin{figure}[tph]
\begin{center}
\input{budgetTF.tex}
\end{center}
\caption{Trends of TF, EDF, and AGCF }\label{FIG:sfB}
\end{figure}

We further analyze the trend of EDF by looking into its two components: EDUPF and EDDNF, and we show the \zb{trends} of these three factors in Fig. \ref{FIG:edfB}.  From the \lgreen{figure,} we can observe that EDDNF does not change as the budget changes but \lred{EDUPF} increases and becomes stable \lgreen{when  a higher budget is available}. That is because the operational cost only includes the fuel cost of thermal units, which means that it can reflect the cost caused by increasing the thermal \lgreen{generators'} output. In other words, the economic dispatch \lgreen{upward} flexibility can be restricted by the system budget when the system experiences a high load and a low renewable energy output. Therefore, the budget of the operational cost has a significant impact on EDUPF, especially when the budget is close to the cost of the economic dispatch only without dynamic AGC. On the contrary,  since there is no cost when generators reduce their output, the budget has little impact on the economic dispatch \lgreen{downward} flexibility.
\begin{figure}[tph]
	\begin{center}
	\input{budgetEDF.tex}
\end{center}
\caption{Trends of EDF, EDUPF, and EDDNF }\label{FIG:edfB}
\end{figure}

In addition, we show the trends of AGCF, AGCUPF, and AGCDNF in Fig. \ref{FIG:agcB}. We can observe that the AGCF is monotonic non-decreasing as \lgreen{the} budget increases.  Similarly, AGCUPF and AGCDNF both increase (not monotonic) as \lgreen{the} budget increases.  Once the scale factor of the budget reaches B18, all flexibility indicators remain as constants.

\begin{figure}[tph]
	\begin{center}
\input{budgetAGCF.tex}
\end{center}
\caption{Trends of AGCF, AGCUPF, and AGCDNF }\label{FIG:agcB}
\end{figure}

\subsection{Impacts of Ramping Capability}
In this subsection, we report simulation results by increasing the ramping capability of generators under three different budget settings (i.e., B0, B14, and B19).  Then, we test different ramping capability settings of all generators in the system to check how flexibility indicators change with ramping capability under a fixed budget.   We report  flexibility indicators (i.e., TF, EDF, and AGCF) in Table \ref{tab:SFUDRC} and show the TF curves in Fig. \ref{FIG:rampTF}.

Fig. \ref{FIG:rampTF} shows that the increasing ramping capability can improve both TF and EDF of the system, but it has \lgreen{a} limited effect on enhancing AGCF. Under both budgets (B0 and B14), AGCF \zb{increases} as the SF increases, while under budget B19, the AGCF remains the same with the change of ramping capability. Therefore, we can observe that the ramping capacity will \lred{influence} the economic dispatch significantly but \lred{influence} the dynamic AGC \lred{minimally}.
\begin{table}[h]
	\huge
	\centering
	\caption{{System Flexibility under Different Ramping Capabilities} }
	\label{tab:SFUDRC}
	\begin{adjustbox}{max width=0.48\textwidth}
		\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
			\hline
			\multirow{2}{*}{SF} & \multicolumn{3}{c|}{B0}   & \multicolumn{3}{c|}{B14}  & \multicolumn{3}{c|}{B19}  \\ \cline{2-10}
			& TF      & EDF     & AGCF  & TF      & EDF     & AGCF  & TF      & EDF     & AGCF  \\ \hline
			1                   & 115.298 & 115.298 & 0.000 & 130.878 & 129.476 & 1.402 & 131.976 & 129.476 & 2.500 \\ \hline
			1.01                & 115.431 & 115.431 & 0.000 & 130.904 & 129.501 & 1.403 & 132.001 & 129.501 & 2.500 \\ \hline
			1.02                & 115.564 & 115.563 & 0.000 & 130.930 & 129.526 & 1.404 & 132.026 & 129.526 & 2.500 \\ \hline
			1.03                & 115.697 & 115.696 & 0.000 & 130.956 & 129.551 & 1.405 & 132.051 & 129.551 & 2.500 \\ \hline
			1.04                & 115.830 & 115.829 & 0.001 & 130.982 & 129.576 & 1.406 & 132.076 & 129.576 & 2.500 \\ \hline
			1.05                & 115.963 & 115.962 & 0.001 & 131.007 & 129.601 & 1.406 & 132.101 & 129.601 & 2.500 \\ \hline
			1.06                & 116.096 & 116.095 & 0.001 & 131.033 & 129.626 & 1.407 & 132.126 & 129.626 & 2.500 \\ \hline
			1.07                & 116.229 & 116.227 & 0.001 & 131.059 & 129.651 & 1.408 & 132.151 & 129.651 & 2.500 \\ \hline
			1.08                & 116.362 & 116.360 & 0.001 & 131.085 & 129.676 & 1.409 & 132.176 & 129.676 & 2.500 \\ \hline
			1.09                & 116.494 & 116.493 & 0.001 & 131.111 & 129.701 & 1.410 & 132.201 & 129.701 & 2.500 \\ \hline
			1.1                 & 116.623 & 116.622 & 0.002 & 131.136 & 129.726 & 1.410 & 132.226 & 129.726 & 2.500 \\ \hline
			1.2                 & 117.906 & 117.902 & 0.003 & 131.394 & 129.976 & 1.418 & 132.476 & 129.976 & 2.500 \\ \hline
			1.3                 & 119.152 & 119.147 & 0.005 & 131.652 & 130.226 & 1.426 & 132.726 & 130.226 & 2.500 \\ \hline
			1.4                 & 120.367 & 120.361 & 0.006 & 131.909 & 130.476 & 1.433 & 132.976 & 130.476 & 2.500 \\ \hline
			1.5                 & 121.573 & 121.566 & 0.006 & 132.166 & 130.726 & 1.440 & 133.226 & 130.726 & 2.500 \\ \hline
			1.6                 & 122.746 & 122.739 & 0.007 & 132.423 & 130.976 & 1.447 & 133.476 & 130.976 & 2.500 \\ \hline
			1.7                 & 123.864 & 123.857 & 0.007 & 132.680 & 131.226 & 1.454 & 133.726 & 131.226 & 2.500 \\ \hline
			1.8                 & 124.962 & 124.955 & 0.008 & 132.936 & 131.476 & 1.460 & 133.976 & 131.476 & 2.500 \\ \hline
			1.9                 & 126.054 & 126.046 & 0.008 & 133.193 & 131.726 & 1.467 & 134.226 & 131.726 & 2.500 \\ \hline
			2                   & 127.121 & 127.113 & 0.008 & 133.449 & 131.976 & 1.473 & 134.476 & 131.976 & 2.500 \\ \hline
			2.5                 & 131.161 & 131.151 & 0.011 & 133.775 & 132.259 & 1.515 & 134.760 & 132.259 & 2.500 \\ \hline
			3                   & 132.315 & 132.259 & 0.056 & 133.811 & 132.259 & 1.552 & 134.760 & 132.259 & 2.500 \\ \hline
			4                   & 132.393 & 132.259 & 0.133 & 133.874 & 132.259 & 1.615 & 134.760 & 132.259 & 2.500 \\ \hline
		\end{tabular}
	\end{adjustbox}
\end{table}

\begin{figure}[tph]
	\begin{center}
		\input{rampTF.tex}
	\end{center}
	\caption{Trends of TF under Different Ramping Capabilities }\label{FIG:rampTF}
\end{figure}

\subsection{Impacts of Transmission Line Capacity}
In this subsection, we study the impact of transmission line capacity on the system flexibility. Similar \lgreen{to} the setup of ramping capability experiments,  we study the system flexibility under three different budgets (i.e., B0, B14, and B19). Then we test different scenarios of all the transmission line capacity in Table \ref{tab:SFUDTLC}  by changing the scale factor.

We can observe that under the budget B0 and B19, AGCF remains unchanged, and under budget B14, it \lgreen{changes} slightly. \zb{This} also indicates that improving the transmission capacity \zb{has a} marginal impact on the AGCF. On the contrary, Fig. \ref{FIG:lineTF} demonstrates that a larger  transmission line capacity can improve TF for a fixed budget until TF reaches its maximum value.
\begin{table}[h]
	\huge
	\centering
	\caption{{System Flexibility under Different Transmission Line Capabilities} }
	\label{tab:SFUDTLC}
	\begin{adjustbox}{max width=0.48\textwidth}

	\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
		\hline
		\multirow{2}{*}{SF} & \multicolumn{3}{c|}{B0}   & \multicolumn{3}{c|}{B14}  & \multicolumn{3}{c|}{B19}  \\ \cline{2-10}
		& TF      & EDF     & AGCF  & TF      & EDF     & AGCF  & TF      & EDF     & AGCF  \\ \hline
		1                   & 115.298 & 115.298 & 0.000 & 130.878 & 129.476 & 1.402 & 131.976 & 129.476 & 2.500 \\ \hline
		1.01                & 115.560 & 115.560 & 0.000 & 131.076 & 129.673 & 1.402 & 132.173 & 129.673 & 2.500 \\ \hline
		1.02                & 115.823 & 115.823 & 0.000 & 131.273 & 129.870 & 1.403 & 132.371 & 129.870 & 2.500 \\ \hline
		1.03                & 116.054 & 116.054 & 0.000 & 131.470 & 130.068 & 1.403 & 132.568 & 130.068 & 2.500 \\ \hline
		1.04                & 116.251 & 116.251 & 0.000 & 131.667 & 130.265 & 1.403 & 132.765 & 130.265 & 2.500 \\ \hline
		1.05                & 116.448 & 116.448 & 0.000 & 131.865 & 130.462 & 1.403 & 132.962 & 130.462 & 2.500 \\ \hline
		1.06                & 116.645 & 116.645 & 0.000 & 132.062 & 130.659 & 1.403 & 133.159 & 130.659 & 2.500 \\ \hline
		1.07                & 116.842 & 116.842 & 0.000 & 132.259 & 130.856 & 1.403 & 133.356 & 130.856 & 2.500 \\ \hline
		1.08                & 117.038 & 117.038 & 0.000 & 132.456 & 131.053 & 1.403 & 133.554 & 131.053 & 2.500 \\ \hline
		1.09                & 117.217 & 117.217 & 0.000 & 132.653 & 131.251 & 1.403 & 133.751 & 131.251 & 2.500 \\ \hline
		1.1                 & 117.392 & 117.392 & 0.000 & 132.851 & 131.448 & 1.403 & 133.948 & 131.448 & 2.500 \\ \hline
		1.2                 & 118.913 & 118.913 & 0.000 & 133.672 & 132.259 & 1.413 & 134.760 & 132.259 & 2.500 \\ \hline
		1.3                 & 119.008 & 119.008 & 0.000 & 133.683 & 132.259 & 1.424 & 134.760 & 132.259 & 2.500 \\ \hline
		1.4                 & 119.008 & 119.008 & 0.000 & 133.687 & 132.259 & 1.428 & 134.760 & 132.259 & 2.500 \\ \hline
		1.5                 & 119.008 & 119.008 & 0.000 & 133.687 & 132.259 & 1.428 & 134.760 & 132.259 & 2.500 \\ \hline
		1.6                 & 119.008 & 119.008 & 0.000 & 133.687 & 132.259 & 1.428 & 134.760 & 132.259 & 2.500 \\ \hline
		1.7                 & 119.008 & 119.008 & 0.000 & 133.687 & 132.259 & 1.428 & 134.760 & 132.259 & 2.500 \\ \hline
		1.8                 & 119.008 & 119.008 & 0.000 & 133.687 & 132.259 & 1.428 & 134.760 & 132.259 & 2.500 \\ \hline
		1.9                 & 119.008 & 119.008 & 0.000 & 133.687 & 132.259 & 1.428 & 134.760 & 132.259 & 2.500 \\ \hline
		2                   & 119.008 & 119.008 & 0.000 & 133.687 & 132.259 & 1.428 & 134.760 & 132.259 & 2.500 \\ \hline
		2.5                 & 119.008 & 119.008 & 0.000 & 133.687 & 132.259 & 1.428 & 134.760 & 132.259 & 2.500 \\ \hline
		3                   & 119.008 & 119.008 & 0.000 & 133.687 & 132.259 & 1.428 & 134.760 & 132.259 & 2.500 \\ \hline
		4                   & 119.008 & 119.008 & 0.000 & 133.687 & 132.259 & 1.428 & 134.760 & 132.259 & 2.500 \\ \hline
	\end{tabular}

	\end{adjustbox}
\end{table}

\begin{figure}[tph]
	\begin{center}
		\input{lineTF.tex}
	\end{center}
	\caption{Trends of TF under Different Transmission Line Capabilities }\label{FIG:lineTF}
\end{figure}

\section{{Conclusions}}\label{sec:conclu}
In this paper, we proposed a real-time flexibility management framework to model economic dispatch with dynamic AGC constraints. The proposed framework can be solved by conducting reformulation and decomposition. We further proposed \lred{seven system flexibility indices} (i.e., TF, EDF, AGCF, EDUPF, EDDNF, AGCUPF, AGCDNF) to reflect the system flexibility and studied how the system factors such that the operational cost budget, ramping capacity and transmission line capacity can impact the system flexibility. We found that the budget of the operational cost can \zb{significantly} contribute to all indicators expect EDDNF, and the improvement of ramping capability can significantly enhance the EDF and the AGCF. In addition, we discover that the transmission capacity only \lgreen{contributes} to EDF.  As the future work, we will study the \zb{power system} flexibility with considering the inertia of the resources in the system.

\bibliographystyle{IEEEtran}
\bibliography{CC_Edge}

\end{document}

config: arxiv


text:
\section{Introduction.}

Examining the behavior of fluid flow in elementary domains is more than a
common numerical test for a new fluid model, it is also a means of
highlighting interesting phenomena inherent in the model.
The direct numerical simulation (DNS) of turbulent
flow at small to moderate Reynolds number has become an important
computational tool in understanding large scale turbulence
motion.  Unfortunately, DNS are still computationally expensive in
turbulent regimes.  Alternate
approaches to brute-force DNS are the Reynolds averaged
Navier-Stokes (RANS) simulations and large eddy simulations (LES).  The
RANS equations (\ref{eq:RANS}) are briefly discussed in Chapter
\ref{chapter:introduction}.
In LES, a spatial averaging operator (filter) is applied to the
Navier-Stokes equations to obtain a new set of equations for the
averaged (filtered) variables.  Due to the computational limit of DNS
for large Reynolds number, LES have become one of the
standard methods in solving for fluid flow.

The behavior of small spatial scales in
turbulent flow is often characterized by statistical isotropy
and homogeneity away from the boundary of the fluid container.
Therefore, the
isotropic LANS-$\alpha$ equations
appear to be an appropriate model in isotropic regimes and have
been studied recently from the numerical point of view.  As a proposed model
for large scale turbulence, Mohseni \emph{et. al.} \cite{msk03}
compare the isotropic LANS-$\alpha$ equations to known results for
DNS and LES methods.  The authors demonstrate the utility of the isotropic
~LANS-$\alpha$ as a sub-grid stress model for three-dimensional isotropic
forced and decaying turbulence.  They perform two sets of
forced isotropic turbulence simulations and compare their results to
DNS and LES results where appropriate.   In the LANS-$\alpha$
simulations, Mohseni \emph{et. al.} conclude that
selecting an appropriate $\alpha$ is a delicate
comprise between the accuracy of the model at the the large scales and
the minimum resolution requirements.  In particular, the
LANS-$\alpha$ equations accurately mimic the behavior of the
Navier-Stokes equations at large spatial scales as long as a minimum
resolution is observed.  The accuracy improves as $\alpha \rightarrow
0$,
but the computation requires higher resolution.  The higher resolution
comes with a price: an unresolved LANS-$\alpha$ computation could result in
the contamination of large scales and a loss of accuracy.

In 1998, Chen \emph{et. al.} \cite{cfh98} studied the mean velocity of
turbulent channel and pipe flows.  They proposed using stationary
solutions to isotropic LANS-$\alpha$ equations as a closure
approximation for the Reynolds-averaged equations.  Since in the
near-wall region the fluctuations are highly anisotropic, their
results were only in good agreement with experimental data away from the
viscous boundary layer.

By restricting the anisotropic LANS-$\alpha$ equations to the
channel, Coutand and Shkoller \cite{cs03a} propose a turbulent channel theory
that models the large scale fluid motion throughout the entire domain.  Unlike
the isotropic equations, the solutions to the anisotropic equations consist
of both a mean velocity field $u$ and a $3 \times 3$ covariance tensor
$F$, defined, in detail, in the next section.  An important property of $F$
is its degeneracy to zero at the wall.   Coutand
and Shkoller \cite{cs03a} show that near the
wall, the degeneracy scales like $d \sqrt{ | \log d |}$, where $d$ is the
normalized distance function to the wall.  The authors
compensate for the degeneracy in the boundary layer by working in
weighted Sobolev spaces.  In this functional framework, they prove the
global-in-time existence and uniqueness of weak solutions to the
anisotropic LANS-$\alpha$ equations.  In particular, they restrict the
anisotropic equations to the channel and make the
assumption that the initial covariance tensor is given in the form
$F_{0}=F(0,x,y,z) = \rho(z) \, \mbox{Id}$.  By assuming the fluid is moving in
one direction, the anisotropic LANS-$\alpha$ equations reduce
to a one-dimensional partial differential equation for the mean velocity.
The authors use the Galerkin method to obtain the existence of weak solutions
to the anisotropic equations.

In this chapter, we use the anisotropic LANS-$\alpha$ equations as a numerical
model for two classical examples of laminar velocity profiles.   Unlike
other approaches to the modeling of tubulence, the Lagrangian averaging
approach allows us to capture the large scale motion of a fluid in laminar
regimes.  The
first example that we consider is the Poiseuille flow in the channel and pipe
domains.  Laminar Poiseuille flow occurs when an incompressible fluid
with no-slip boundary conditions is driven by a constant upstream pressure
gradient, yielding a symmetric parabolic stream-wise profile.   We assume
the velocity $u$ satisfies the steady Navier-Stokes equations and we
numerically solve
the anisotropic equations for the initial covariance tensor $F_{0}$.  In
particular, we calculate the matrix $F$ such that pair $(u,F)$ is a
solution to the
anisotropic LANS-$\alpha$ equations.  The degeneracy rate of our
numerical solution $F$ near the wall of the channel is in good agreement with the logarithmic decay rate
given by Coutand and Shkoller \cite{cs03a}.  In addition, we show that in the
boundary layer, $||F(t, \cdot)||_{L^{\infty}(\Omega)}$ does not remain
bounded as $t \rightarrow \infty$, answering, at least numerically, a question
stated in \cite{cs03a}.

In the last section of the chapter, we study shear flow solutions to
the anisotropic LANS-$\alpha$ equations.  Shear flow occurs in the channel
when one side of the boundary is moving while the other side remains fixed.
The derivation by Marsden and Shkoller \cite{ms03} relied
upon no-slip boundary conditions for the mean velocity field $u$.  In order to
find a solution to the inhomogeneous problem, we use the classical
technique of introducing a new vector field
that vanishes on the boundary and that solves the anisotropic equations with
additional forcing terms.  For the domains and the initial data that submit a
unique solution to the new formulation, we are able to solve for the
unknown mean velocity.  In Section \ref{inhomog}, we show that shear flow
velocity solutions to the anisotropic LANS-$\alpha$ equations exist if the
initial covariance tensor is not required to be positive.  This turns out
to be a natural assumption and we compute the shear flow solutions directly.

\subsection{The covariance tensor $F(t,x).$}

The matrix $F$ is defined by Marsden and Shkoller
\cite{ms03}
as the ensemble average of the tensor product of the Lagrangian
fluctuation vector.  Specifically, if $\eta$ and $\eta^{\varepsilon}$
are the particle trajectories of the fluid velocity $u$ and averaged
fluid velocity $u^{\varepsilon}$, respectively, then the Lagrangian
fluctuation vector
$\xi^{\varepsilon} := \eta^{\varepsilon} \circ \eta^{-1}$ (see Figure
\ref{fig:xi}) and
\[ F := \displaystyle{ \left<
\frac{d}{d\varepsilon}\Big|_{\varepsilon=0} \xi^{\varepsilon} \otimes
\frac{d}{d\varepsilon}\Big|_{\varepsilon=0} \xi^{\varepsilon} \right> }\]
where $<\cdot>$ denotes the average over all possible solutions
$u^{\varepsilon}$ (see \cite{ms03}).
\begin{figure}[tbp]
    \centering
    \resizebox{10cm}{!}{\includegraphics{xi.eps}}
    \caption{The construction of $F$.}
    \label{fig:xi}
\end{figure}
Because these fluctuations are necessarily zero along the boundary,
\[ F(t,x) = 0 \,\, \mbox{for} \,\, t \geq 0 \]
on the boundary of the domain.  As noted in \cite{cs03a}, it is
unknown \emph{a priori} whether
$||F(t, \cdot)||_{L^{\infty}}$ remains bounded for all time.
In next section, we show numerically that $||F(t, \cdot)||_{L^{\infty}}$ is
strictly increasing in time in the viscous boundary layer of the channel and
pipe under a steady mean flow.

\section{Homogeneous boundary conditions.}\label{homog}

Laminar Poiseuille flow occurs when an incompressible fluid in a straight
channel, or pipe, is driven by a constant upstream pressure gradient, yielding a
symmetric parabolic stream-wise velocity profile.  In this section,
we begin with the steady Poiseuille flow associated with the Navier-Stokes
equations.  We construct an initial covariance tensor such that the
steady flow solution and the
covariance tensor solve the anisotropic equations in the channel.  First,
we reduce the anisotropic
LANS-$\alpha$ equations to the channel and pipe under the assumption that the
initial covariance tensor is a multiple of the identity matrix.

\subsection{Given the initial covariance matrix.}

\noindent {\bf Channel.}   The three-dimensional channel is given by
$\Omega = \mathbb{R}^2 \times
[-h,h]$ with coordinates ${\bf x}=(x,y,z)$.  We shall assume that the velocity
 vector is of the form
\[ {\bf u}(t,x,y,z) = (u(t,z),0,0), \]
and that the initial covariance matrix is written as
\[ F_{0}=F(0,x,y,z) = \rho(z) \mbox{Id}, \]
such that $\rho(\pm h) = 0$.

Suppose $\rho(z)$ is given.  As discussed in \cite{cs03a}, when no-slip
boundary conditions are prescribed for the mean velocity, the
nonlinear term vanishes and equation (\ref{aniso}) reduces to the
following system
\begin{equation}\label{aniso_channel}
    \begin{array}{c}
\partial_t (u-\alpha^2(\rho u')')-\nu((\rho u')' -
\alpha^2 (\rho(\rho u')'')') = - c, \\
u(t,\pm h) = 0, \hspace{.2in} u(0,z)=u_0(z),
\end{array}
\end{equation}
where $c :=\partial_x p$ is constant since $\partial_y p= \partial_z p
=0$.  To address the degeneracy at the boundary, Coutand and Shkoller
define a weighted Sobolev to set as the functional framework.  This
space is based on the following function:  For $\epsilon < \delta << 1$
define $\rho_{CS}(z) \in C^{\infty} [-h,h]$ to be the positive function
\begin{equation} \label{p_choosen}
    \begin{array}{c}
\rho_{CS}(z) = \left\{
\begin{array}{cc}
d \sqrt{|\log d|}, & 0 \leq d(z) \leq \epsilon, \\
1, & d(z) \geq \delta,
\end{array}
\right.  \\
d \ \, \mbox{= normalized distance function to the boundary}.
\end{array}
\end{equation}
Coutand and Shkoller prove the existence of a unique global
weak solution to equation (\ref{aniso_channel}) in a weighted
Sobolev space.  In addition, they show that the covariance tensor $F$ must
degenerate like $\rho_{CS}$ in the viscous boundary layer.  In the next section,
we verify that our numerical solution $\rho$ matches $\rho_{CS}$ near
the walls of the channel.

\noindent {\bf Pipe.}  The three-dimensional
pipe is given by $\Omega = \mathbb{R}^3$ with
coordinates ${\bf x} = (x,y, z)$ such that $0\leq \sqrt{x^2+y^2}\leq a$.
We shall assume that the velocity vector is radially symmetric of the form
\begin{equation}\label{u_pipe}
{\bf u}(t,x,y,z) = (0,0,u(t,r)),
\end{equation}
where $r:=\sqrt{x^2+ y^2}$ and that the initial covariance matrix
\[ F_{0}=F(0,x,y,z) = \rho(r) \mbox{Id}.\]
We let ${\bf \eta}$ denote the Lagrangian flow of ${\bf u}$ satisfying the
initial value problem (\ref{eq:eta}).
Then definition (\ref{u_pipe}) implies that
\[ \displaystyle{ {\bf \eta}(t,{\bf x}) = \left(x,y, z+\int_0^t u(s,r) \, ds
\right).} \]
Therefore
\begin{equation} \label{F_pipe}
\displaystyle{ F(t,{\bf x}) := D{\bf \eta}(t,{\bf x}) \cdot F_0({\bf x}) \cdot
D{\bf \eta}(t,{\bf x})^T = \rho(r) \left[
\begin{array}{ccc}
1 & 0 & \frac{x}{r}\mathcal{U} \\
0 & 1 & \frac{y}{r}\mathcal{U} \\
\frac{x}{r}\mathcal{U} & \frac{y}{r}\mathcal{U} & 1+\mathcal{U}^2
\end{array}
\right],}
\end{equation}
where $\displaystyle{\mathcal{U}(t,{\bf x}):=\int_0^t \partial_r u(s,r)\,ds}$.
Together with definition (\ref{u_pipe}), $C{\bf u} := \operatorname{div}(\nabla u \cdot
F)$ reduces to
\[ C{\bf u} = (0,0,\frac{1}{r}\left(r \rho u'\right)') \]
with $u'(t,r):=\partial_r u(t,r)$.  The nonlinear term in the
anisotropic LANS-$\alpha$ vanishes and the equations reduce to the following
system
\begin{equation}\label{aniso_pipe}
    \begin{array}{c}
\displaystyle{\partial_t (u- \alpha^2\frac{1}{r}(r \rho u')')-
\frac{\nu}{r}\left( (r\rho u')'-\alpha^2  (r\rho
(\frac{1}{r}(r\rho u')')')'\right) = - c,} \\
u(t,a) = 0, \hspace{.2in} u(0,r)=u_0(r),
\end{array}
\end{equation}
where $c:=\partial_z p $ is constant since $\partial_x p= \partial_y
p =0$.

\subsection{Given the mean velocity profile.}\label{sec:num_F}

In this section, we provide the numerical
results for steady channel and pipe flow assuming no-slip boundary conditions.
We start with the steady solution $u$ to the Navier-Stokes equations and proceed
to find a $\rho$, and therefore a $F$, such that $(u,F)$ solve the anisotropic
LANS-$\alpha$ equations.

\noindent {\bf Channel.}  We begin by assuming the velocity
field is the steady flow given as
${\bf u}(t,x,y,z) = (u(z),0,0)$ that satisfies the no-slip boundary condition
$u(\pm h) = 0$.  For the Navier-Stokes equations, the classical Poiseuille
flow is given by
\begin{eqnarray}\label{poiseuille}
u(z) &=& \beta (h^2 - z^2) \\
p(x) &=& p_0 - 2\nu \beta x, \nonumber
\end{eqnarray}
where $\beta, p_0 >0$.  In order for $(u,\rho)$ to
solve equation (\ref{aniso_channel}) with the velocity $u(z)$ given
by equation (\ref{poiseuille}), the function $\rho(z)$ must solve the
ordinary differential equation
\begin{equation}\label{channel_rho}
    \begin{array}{c}
(z\rho)'-\alpha^2 (\rho (z \rho)'')' =1, \\
\rho(\pm h) = 0 .
\end{array}
\end{equation}
Solving equation (\ref{channel_rho}) numerically leads us to the following
Proposition.

\begin{proposition} Let the function $u$ be the steady solution
    (\ref{poiseuille}) to the incompressible Navier-Stokes equations in
    the smooth three-dimensional
    channel.  The pair $(u,\rho)$ solves the channel anisotropic
LANS-$\alpha$ equation (\ref{aniso_channel}) when the function $\rho(z)$
solves equation (\ref{channel_rho}).  A numerical solution is given in
Figure \ref{fig:rho} when $\alpha = 0.1$ and $h = \beta = 1$.  In addition,
the covariance tensor $F$ is defined as
\begin{equation*}
F(t,z) := \rho (z) D\eta(t,{\bf x}) D\eta(t,{\bf x})^T = \rho(z) \left[
\begin{array}{ccc}
1+4\beta^2 t^2 z^2 & 0 & -2\beta tz \\
0 & 1 & 0 \\
-2 \beta tz & 0 & 1
\end{array}
\right].
\end{equation*}
\end{proposition}

\begin{remark} We assume that $\rho(z)$ is positive and
    shares the same symmetry across the channel as the velocity $u(z)$.  In
    addition, we assign in the center of the channel the values
    $\rho(0) = A$ and $\rho^{\prime}(0)=0$.  By a numerical shooting method, we
    find a value for $A$ that forces the zero boundary condition $\rho(1)=0$.
    This solution is plotted in Figure \ref{fig:rho}.
\end{remark}

\begin{remark} The numerical solution $\rho$ decays at the same rate as the
    function $\rho_{CS}$ near the boundary.  In Figure \ref{fig:CS_rho},
    we plot both functions near the wall of the channel.
\end{remark}

To illustrate the dynamics of $F$ in the channel, we calculate
its eigenvalues and conclude the following Proposition.

\begin{proposition} \label{F_prop}
    The eigenvalues of $F(t,z)$ are
\begin{eqnarray*}
\lambda_1(t,z)&=& \rho(z), \\
\lambda_{2,3}(t,z) &=& \frac{1}{2} \rho(z) \left( 2+\mathcal{U}^2 \pm
|\mathcal{U}|\sqrt{\mathcal{U}^2+4} \right) = \rho(z) \left(1+2\beta^{2}
t^2z^2 \pm 2\beta t|z| \sqrt{1+\beta^{2} t^2z^2} \right).
\end{eqnarray*}
The eigenvalues are plotted in Figure \ref{fig:ev_of_F} for the case of
$\alpha = 0.1$, $\beta =1$, and the function $\rho$ given as the solution
to equation (\ref{channel_rho}).  Furthermore, $||F(t,\cdot)||_{L^\infty(\Omega)} =
O(t^2)$ and although our velocity is steady throughout the channel, the
covariance tensor
\[ ||F(t,\cdot)||_{L^{\infty}(\Omega)} \rightarrow \infty \hspace{.2in}
\mbox{as} \hspace{.2in} t \rightarrow \infty \]
near the boundary in the direction of the flow.  This is illustrated
in Figure \ref{fig:F_max}.
\end{proposition}

\newpage

\begin{figure}[!ht]
     \centering
     \subfigure[The solution $\rho(z)$ of equation (\ref{channel_rho}). ]{
	  \label{fig:rho}
	  \includegraphics[width=.45\textwidth]{channel1.eps}}
    \hspace{.2in}
     \subfigure[The solution $\rho(z)$ and the function $\rho(z)$ used in
     \cite{ms03} near the boundary.]{
	  \label{fig:CS_rho}
	  \includegraphics[width=.45\textwidth]{channel2.eps}} \\
     \vspace{.2in}

     \subfigure[The eigenvalues of $F(t,z)$ at $t=2$.]{
	   \label{fig:ev_of_F}
	   \includegraphics[width=.45\textwidth]
		{channel3.eps}}
    \hspace{.2in}
     \subfigure[The evolution of $||F||_{L^\infty}$.]{
	   \label{fig:F_max}
	  \includegraphics[width=.45\textwidth]{channel4.eps}}
     \caption{For $h=\beta =1$ and $\alpha=0.1$. Because of symmetry,
     the graphs (a), (c), and (d) are on $[0,1]$ only.}
     \label{fig:channel}
\end{figure}

\newpage

\noindent {\bf Proof.} The eigenvalues are straightforward to compute.  Since
the largest eigenvalue
$\lambda_2(\cdot, x) = O(t^2)$ and $\lambda_2 \leq ||F(t,z)||_{L^\infty}$,
then $||F(t, \cdot)||_{L^\infty} = O(t^2)$.  The vector
\begin{equation*}
v(t,z) = C \left[
\begin{array}{c}
tz + \sqrt{1+t^2 z^2} \\
0 \\
-1
\end{array} \right]
\end{equation*}
solves the eigenvalue equation $Fv=\lambda_2v$.  Therefore, as
$t \rightarrow \infty$,
the norm of $F$ increases to infinity in the direction of the flow
($x$-direction).
\qed

\noindent {\bf Pipe.}  The format and conclusions of this subsection
follow closely the results in the channel.  We begin by assuming the velocity
field is a steady velocity field given as
${\bf u}(t,{\bf x}) = (0,0,u(r))$ that satisfies the no-slip boundary
condition $u(\pm a) = 0$.  For the Navier-Stokes equations, the classical
pipe flow is given by
\begin{equation}\label{pipe}
    \begin{array}{c}
u(r) = \beta (a^2 - r^2) \\
p(r) = p_0 - 4\nu \beta r,
\end{array}
\end{equation}
where $\beta, p_0 >0$.  In order for $(u,\rho)$ to
solve equation (\ref{aniso_pipe}) with the velocity $u(r)$ given
by equation (\ref{pipe}), the function $\rho(r)$ must solve the
ordinary differential equation

\begin{equation}\label{eq:pipe_rho}
   \begin{array}{c}
       r(r\rho'+2\rho) - \alpha^2(r\rho (r\rho''+3\rho'))' = 2r \\
       \rho(a) = 0.
\end{array}
\end{equation}
As in the channel domains, we arrive at the following Proposition.

\begin{proposition} Let the function $u$ be the steady solution (\ref{pipe})
to the incompressible Navier-Stokes equations in the smooth three-dimensional
   pipe.  The pair $(u,\rho)$ solves the pipe anisotropic
LANS-$\alpha$ equation (\ref{aniso_pipe}) when the function $\rho(r)$ solves
equation (\ref{eq:pipe_rho}).  A numerical solution is given in Figure
\ref{fig:pipe_rho} when $\alpha = 0.1$, $\beta =1$, and the radius of the pipe is $a = 1$.
In addition, the covariance tensor $F$ is given by (\ref{F_pipe}) and
has the same eigenvalues as the channel covariance tensor.
\end{proposition}

\begin{figure}[thp]
    \centering
    \resizebox{8cm}{!}{\includegraphics{pipe_rho.eps}}
    \caption{The solution $\rho(r)$ of equation (\ref{eq:pipe_rho}.)}
    \protect\label{fig:pipe_rho}
\end{figure}

\subsection{Conclusions.}

As expected, the anisotropy of the fluid is of fundamental importance
in bounded domains.  We have demonstrated that the dynamics of $F$ are
just as important in laminar regimes as they are in turbulent regimes.
Supposing that the covariance tensor is a constant multiple of
the identity matrix, while accurate for very short time and in the
center of the channel, becomes an inaccurate assumption in the viscous
boundary layer.  In both the channel and the pipe geometry, the evolution of
$F$, illustrated by Figure \ref{fig:ev_of_F}, is non-decreasing in
time and achieves its greatest value in this boundary layer.   As time
increases, the location of the maximum value of
$||F||_{L^{\infty}}$ in the channel approaches the limit $z=0.89$.
Therefore an accurate model for fluid motion in the entire channel or
pipe should be founded on the anisotropic model, rather than the isotropic
version, in this boundary region.

The logarithmic degeneracy rate of the covariance tensor in the channel
agrees well with the decay rate of the function computed in \cite{cs03a}.
Since the dynamics of fluid motion in the channel and pipe are similar, it
is not surprising that the covariance tensor for the channel and pipe are the
same.  The
next step in studying the numerical properties of the anisotropic model is
to solve for the covariance tensor when the velocity is given as a
time-dependent solution to the Navier-Stokes problem.  This is a necessary
step in understanding the dynamics of the covariance tensor $F$.

\section{Inhomogeneous boundary conditions.} \label{inhomog}

In many physical models, ranging from the study of blood flow to the
modeling of earthquakes, at least one of the boundary components of
the fluid container are in motion.  In these cases,
the mean velocity of the fluid at the wall will no longer be zero.  In this
section, we study the anisotropic LANS-$\alpha$ equations with
inhomogeneous boundary data.  As a specific example, we consider the mean
fluid motion in a channel when one of the boundary walls is not fixed.
This motion is called shear flow.

\subsection{The LANS-$\alpha$ equations with inhomogeneous boundary
data.}

Let $\Omega$ be a bounded fluid container in $\mathbb{R}^{n}$ with
boundary $\partial \Omega$ and suppose that the mean velocity field
$u=g \neq 0$ on $\partial \Omega$.  To find a solution pair $(u,F)$
that solves the anisotropic LANS-$\alpha$ equations with
inhomogeneous boundary conditions, we choose a divergence-free vector
field $v(t, \cdot) \in C^{\infty}(\Omega)$ such that $v=g$ on the
boundary.  With $v$ given, we define a new vector field $w=u-v$.  The
vector field $w$ is zero on the boundary and solves the anisotropic
LANS-$\alpha$ equations with $u$ replaced with $w+v$.  Namely, we now
search for a solution to the following system of partial differential
equations
\begin{equation}\label{aniso_inhom}
    \begin{array}{c}
(1-\alpha^2 C) \left( \partial_t w - \nu \mathbb{P} Cw +\operatorname{div}
(v \otimes w + w \otimes v) \right) = - \operatorname{grad} p - (1-\alpha^2
C) \left( \partial_t v - \nu \mathbb{P} Cv \right), \\
\operatorname{div} w(t,x) = 0, \hspace{.2in} w=0 \,\, \mbox{on} \,\,
\partial \Omega,  \\
\partial_t F + \nabla F \cdot w - \left( \nabla w \cdot F +
[\nabla w \cdot F]^T \right) = \nabla F \cdot v + \nabla v \cdot F +
[\nabla v \cdot F]^T,  \\
w(0,x) = u_0(x)-v(0,x),  F(0,x)= F_0(x).
\end{array}
\end{equation}

In the three-dimensional torus $\mathbb{T}^{3}$, the trivial existence
and uniqueness of a solution to (\ref{aniso_inhom}) follows from the
smoothness of $v$, Theorem 1 in \cite{ms03}, and Theorem
\ref{stokes_wp} in Chapter \ref{chapter:spwp}.  This result is stated
as the following proposition.
\begin{proposition}
    For $s>7/2$, and $w_0 \in
    H^s_{\operatorname{div}}
    (\mathbb{T}^3)$, $F_0 \in [H^s_{\operatorname{per}}
    (\mathbb{T}^3)]^{3 \times 3}$, with $F_0>0$, there exists a unique solution
    $(w, F)$ with $w \in C^0([0,T]; H^s_{\operatorname{div}}) \cap L^2 (0,T;
    H^{s+1}_{\operatorname{div}})$ and $F \in C^0([0,T];
    [H^s_{\operatorname{per}}]^{3 \times 3})$ to equations
    (\ref{aniso_inhom}), where $T$ depends on the initial data.
\end{proposition}

In arbitrary bounded domains, it is unknown whether solutions
exist to the anisotropic equations.  As we demonstrate in the next
section, under certain limiting conditions we can find a solution to
equation (\ref{aniso_inhom}).

\subsection{Given the initial covariance matrix.}

To study shear flow velocity solutions to the anisotropic
LANS-$\alpha$ equations, we need to restrict the full equations to the
three-dimensional channel by making a number of
limiting assumptions.  We assume that the initial covariance tensor $F_0 =
F(0,x,y,z)=\rho(z) \mbox{Id}$ is given and the mean velocity has the
form ${\bf u} = (u(t,z),0,0)$ for $z\in[-h,h]$, where $u \neq 0$ on the
boundary.

As was done in the general case above, we choose a vector
field $v(t,\cdot) \in C^\infty [-h,h]$ such that $v(t,\pm h) = u(t,\pm
h)$ for all $t\geq 0$.  The vector field $w=u-v$ then solves the
following one-dimensional problem
\begin{equation}\label{aniso_channel_inhom}
\begin{array}{c}
    \mathcal{L}^\alpha w = -c +f, \\
w(\pm h) = 0,  w(0,z)=u_0(z)-v(0,z)
\end{array}
\end{equation}
where $\mathcal{L}^\alpha$ is the linear operator defined as
\begin{equation*}
\mathcal{L}^\alpha u := \partial_t(u-\alpha^2(\rho u')')-\nu((\rho u')' -
\alpha^2 (\rho (\rho u')'')'),
\end{equation*}
and $f:= -\mathcal{L}^\alpha v$.  Then ${\bf u}(t,x,y,z) :=
(w(t,z)+v(t,z), 0,0)$ solves equation
(\ref{aniso_channel}) with inhomogeneous boundary conditions.

\subsection{Given the mean velocity field.}

Suppose that we are given the steady velocity vector
field ${\bf u} = (u(z),0,0)$ with $u(z) \neq 0$ on the boundary.
To find a solution to the anisotropic model, we need to find a
function $\rho(z)$ such that
\begin{equation*}
    \nu ((\rho u')'-\alpha^2(\rho (\rho u')'')') = -c.
\end{equation*}
But since the vector field $u$ doesn't satisfy the no-slip boundary
conditions, we do not know the correct boundary conditions for $\rho$.
Rather, we choose $v\in C^{\infty}[-h,h]$ with $v=u$ on boundary.
Then $w=u-v$ and the pair $(w,\rho)$ solve
\begin{equation}\label{channel_inhom}
\mathcal{L}^\alpha w = -c +f,
\end{equation}
where $f:= - \mathcal{L}^\alpha v$.  Since $w=0$ on the boundary we also have
the boundary condition $\rho(\pm h) =0$.
Since $w$ is given we may find
$\rho$ satisfying equation (\ref{channel_inhom}) with zero boundary
conditions.  However, unlike the case when $u$ satisfies no-slip boundary
conditions, the forcing in equation (\ref{channel_inhom}) removes any
a priori statement about the positivity of $\rho$.

\noindent {\bf Shear flow solution.}  As an example, we show the existence of
a shear flow velocity solution to the anisotropic LANS-$\alpha$ equations.  The
steady shear flow velocity solution
to the Navier-Stokes equations with the boundary conditions $u(-h)=0$
and~$u(h)=U>0$ is
\begin{equation}\label{eq:u_shear}
    \displaystyle{u(z) = \frac{U}{2h}(z+h)}
\end{equation}
with the pressure function $p \equiv 0$.  The function
$$\displaystyle{w(z) = \frac{U}{4h^2}(h-z)(h+z)}$$ is zero on the boundary and
solves (\ref{channel_inhom}) with $\displaystyle{f = - \mathcal{L}^\alpha
\left( \frac{U}{4h^2}(z+h)^2 \right)}$ when $\rho(z)$ is a solution of
\begin{eqnarray*}
\rho' -\alpha^2 (\rho \rho'')' &=& 0, \\
\rho(\pm h)&=&0.
\end{eqnarray*}
Namely,
\[ \rho(z) = \left\{
\begin{array}{ll}
\frac{1}{2\alpha^2}(z-a)(z-b) & \mbox{for}\,\, -h \leq a \leq z \leq b \leq h\\
0 & \mbox{otherwise}
\end{array} \right. . \]
If we want $\rho \in C^1[-h,h]$ then $a=-h$, $b=h$.  We conclude that
the pair $(u,\rho)$ with $u$ define by equation (\ref{eq:u_shear})
and $\displaystyle{\rho(z) = \frac{1}{2\alpha^{2}} (z-h)(z+h)}$ is a shear flow
solution which solves equation (\ref{aniso_channel}) with inhomogeneous
boundary conditions.


TITLE: Solve Non-lineaer Differential Equation (Step-by-Step)
QUESTION [0 upvotes]: Being able to solve Linear and Non-liner differentional equations is an important part of mathematics. I am trying to work out the exact solution for this non-linear differential equation and require some help from the community (Please could you provide the full step-by-step-solution), I want to undertstand the mathematics and be able to apply this do other questions.
$$\dfrac {du}{dx}+ u^2 = 0  \tag 1 $$
Boundary conditions $0<x<1; u(0)=1$
What is the exact solution?

REPLY [1 votes]: This is an answer to your question in the comments. Please, keep in mind for the future, provide any attempts you have made, any progress, everything you have tried to solve the problem. Where did you stuck, what was the difficulty? MSE is not a place to simply solve your homework for you ...
Regarding your question about the solution of $u^\prime + u = 0$. Apply the same logic as @Rem explained and also here, think about the problem first:
You want to find a function such that adding its own derivative would results in zero, i.e., they would cancel each other out. So, we need two identical terms with an opposite sign. How about something like -1 + 1 = 0? It seems like we are looking for a function which is very similar to itself when differentiating. Do you have a clue? If not, you can solve this ODE by separation of variables:
First, write your ODE in Leibnitz notation
$$
\frac{\mathrm{d}u}{\mathrm{dx}} + u = 0 \\
\frac{\mathrm{d}u}{\mathrm{dx}} = -u.
$$
We can already see that we are looking for a function which after differentiation stays the same, only changes its sign. Still no clue yet? Change the variables such that the terms containing the variable $u$ are on one side, and the terms containing the dependent variable $x$ on the other side:
$$
\frac{1}{\mathrm{d}u} = -1 \mathrm{d}x.
$$
Now you want to integrate both sides over $x$, so we change the dependent variable from $x$ to $s$ to avoid confusion and write
$$
\int_0^x\frac{1}{u}\mathrm{d}u = \int_0^x-1\mathrm{d}s.
$$
This is easy, think about what function differentiated results in the fraction you see on the l.h.s? Exactly, it's the logarithmic function. If $u = \ln(u)$, then $u^\prime = \frac{1}{u}$, hence we have
$$
\ln(u) + c_1 = -x + c_2.
$$
Bring the coefficient $c_1$ to the right side, apply the exponential function on both sides, et voilà, there's your solution:
$$
%\mathrm{e}^{\ln(u)} = \mathrm{e}^{-x} \cdot \underbrace{\mathrm{e}^{c_1 - c_2}}_{= c_3}\\[1em]
u(x) = \mathrm{e}^{-x} + c_3.
$$
If you want you can check your solution, just take the derivative of this function and add it to the function $u(x)$: $\mathrm{-e}^{-x}+\mathrm{e}^{-x} = 0$.
I forgot, you are given the boundary condition $u(0) = 1$. Apply this to your function and solve for the constant $c_3$:
$$
1 = \mathrm{e}^{-0} + c_3 \quad \Longleftrightarrow \quad c_3 = 0.
$$


TITLE: $\mathbb{N} \in \mathbb{N}$?
QUESTION [1 upvotes]: My professor proved this:
$\forall n \in \mathbb{N}: n \not\in n$
I proved that: $ n \bigcup \{n\} \not\subset n$
pf.
$Suppose,  n \bigcup \{n\} \subset n$
$\implies  (n \bigcup \{n\})\setminus n  \subset n\setminus n$
$\implies  n\setminus n \bigcup \{n\}\setminus n  \subset n\setminus n$
$\implies \emptyset \bigcup \{n\}\setminus n  \subset \emptyset$
$\implies\{n\}\setminus n  \subset \emptyset$
$Suppose, \{n\}\setminus n = \emptyset$
$\implies n \in n$
$ Since, n \in \mathbb{N}, n \not\in n$
By contradiction $\{n\}\setminus n \not= \emptyset$
Therefore, $\{n\}\setminus n = \{n\} \implies\{n\}  \subset \emptyset$
$\implies n \in\emptyset \implies false$
By contradiction $ n \bigcup \{n\} \not\subset n$∎
But, is there any set $X: X\in X$ is true. In particular are the natural numbers in the natural numbers given that they're infinite.

REPLY [3 votes]: But, is there any set $X:X\in X$ is true.

No (at least not in normal set theories, others exist). This is important, and essentially the reason why axiomatic set theory exists: allowing things like this makes incredibly weird things happen. This is called well-foundedness (and it's a direct consequence of the Axiom of Regularity in ZF, the most common set theory).

In particular are the natural numbers in the natural numbers given that they're infinite.

No. If you want to define natural numbers as sets, the usual definition is that $0 = \emptyset$, and $n = \{n - 1\} \cup (n - 1)$. In particular, all such sets are finite, so $\mathbb{N}\not\in\mathbb{N}$.

REPLY [2 votes]: There is nothing inherently wrong with sets $X$ satisfying $X\in X$. They can be reasoned about just like any other set (which is all we can ask for, I guess).
However, one of the axioms in ZF set theory (which is the most commonly used theory, especially if you count it together with ZFC) is called the axiom of regularity. It implies that there are no such sets.
So there is nothing a priori problematic with such sets. They just never occur in the standard theory.


text:
\subsection{Proof of Proposition \ref{prop:trade-off-between-wlc-and-rademacher-complexity}}\label{sec:proof-of-trade-off-between-wlc-and-rademacher-complexity}

\begin{proof}
We denote $p_0 = \frac{MU}{(M+U)^2}$.
Let $\sigma_1, \ldots, \sigma_N$ be the i.i.d. transductive Rademacher variable of parameter $p_0$.
Since $\{-1, 0, 1\}^N\subset \mathcal{V}_{g}$, for any realization of $\bm{\sigma} = (\sigma_1, \ldots, \sigma_N)$, we have $\bm{\sigma} \in \mathcal{V}_{g}$.
By the assumption, there exists $Z_{\bm{\sigma}}\in\mathcal{V}$ such that
\begin{equation*}
    \|Z_{\bm{\sigma}} - \alpha\bm{\sigma}\|_2 \leq \beta \|\bm{\sigma}\|_2.
\end{equation*}
Set $C_f := (2\alpha)^{-1}$ and $C_\mathcal{L} := \frac{\alpha^2 - \beta^2}{2\alpha}$.
Similarly to the proof of Theorem \ref{thm:optimization}, we have
\begin{align*}
    & \|Z_{\bm{\sigma}} - \alpha \bm{\sigma}\|_2 \leq \beta \|\bm{\sigma}\|_2\\
    \Longleftrightarrow \ &
        \|Z_{\bm{\sigma}}\|_{2}^2
            - 2\alpha \langle Z_{\bm{\sigma}}, \bm{\sigma}\rangle
            + \alpha^2 \|\bm{\sigma}\|^2_2
        \leq \beta^2 \|\bm{\sigma}\|^2_2\\
    \Longleftrightarrow \ &
         C_f \|Z_{\bm{\sigma}}\|^2_{2} + C_{\mathcal{L}}\|\bm{\sigma}\|^2_{2}
        \leq \langle Z_{\bm{\sigma}}, \bm{\sigma}\rangle.
\end{align*}
Therefore, we have
\begin{align*}
    Q^{-1} \mathfrak{R}(\mathcal{V})
    &= \E_{\bm{\sigma}} \left[ \sup_{Z\in \mathcal{F}} \langle \bm{\sigma}, Z \rangle \right] \\
    &\geq \E_{\bm{\sigma}} \left[ \langle \bm{\sigma}, Z_{\bm{\sigma}} \rangle \right] \\
    &\geq \E_{\bm{\sigma}} \left[ C_{\mathcal{L}} \|\bm{\sigma}\|^2_2 + C_f \|f_{\bm{\sigma}}\|^2_2 \right]\\
    &\geq \E_{\bm{\sigma}} \left[ C_{\mathcal{L}} \|\bm{\sigma}\|^2_2 \right]\\
    &= C_{\mathcal{L}} \cdot 2Np_0.
\end{align*}
In the last equality, we used \eqref{eq:expectation-of-rademacher-variable}.
Therefore, we have $\mathfrak{R}(\mathcal{F})\geq 2C_{\mathcal{L}}Np_0Q = \frac{\alpha^2 - \beta^2}{\alpha}$.
\end{proof}


TITLE: The function defined by $f(x)=\sin\pi x$ has zeros at every integer. Show that when $-1<a<0$ and $2<b<3$ , the bisection method converges to
QUESTION [0 upvotes]: The function defined by $f(x)=\sin\pi x$ has zeros at every integer. Show that when $-1<a<0$ and $2<b<3$ , the bisection method converges to
(a) $0$, if $a+b<2$
(b) $2$, if $a+b>2$
(c) $1$, if $a+b=2$
$a$ denotes lower limit &  $b$ denotes upper limit of the interval $[a,b]$
How to consider the intervals fulfilling the conditions?
I have thought by the following way:
(a) $1<a+b<2$ . but if i put $f(1)=\sin\pi=0$ and $f(2)=\sin 2\pi=0$
it is not providing the different signs.
I haven't idea to solve the questions.

REPLY [0 votes]: one can notice that $a_0$ lies in the blue region while $b_0$ hence we can establish that : $ 1\leq\frac{a_0 + b_0}{2}\leq\frac{3}{2}$ ; it is clear that after the first iteration if $a_0 + b_0 = 2$ then we get a point $c$ with $f(c) = 0$ and iteration will stop and converge to $1$.
if we say that $a+b < 2$ then we will have the point $c$ for the next iteration to the left of $1$ and the iteration will clearly keep on choosing points $c_n < 1$ and so the iteration converges to $0$.
a similar argument can be checked out for the case when $a+b > 2$ in which iteration converges to $2$ as we always stay on the right of point $(1, 0)$


TITLE: Homomorphism of regular languages
QUESTION [0 upvotes]: I have the homomorphism $h$ of monoid $ \lbrace a,b \rbrace ^* $ which is defined as $h(a) = b , h(b) = aba $.
I know how this is working for example I have in book question is it true that
$h(b^+) \subseteq (a+ba)^* $
and I know that is right answer becouse $h(b^+) = (aba)^+ \subseteq (a+ba)^* $.
But here is another question: is it true that
$h^3(b^+) = (h(b)bh(b)bh(b))^+$?
The answer is yes and I don't know why. I would be very thankful if somebody could tell me how it was done. I don't know why we have here letter $b$ outside the bracket. I have tried some of ways of doing this example but I can't get into this answer $(h(b)bh(b)bh(b))^+$.

REPLY [2 votes]: $h(b^+)=(aba)^+$
$h^2(b^+)=h(h(b^+))=h((aba)^+)=(h(a)h(b)h(a))^+=(bh(b)b)^+$
Therefore $h^3(b^+)=h(h^2(b^+))=h((bh(b)b)^+)=(h(b)h^2(b)h(b))^+  =(h(b)bh(b)bh(b))^+$


TITLE: Show that $a(XY)=X(aY)$ where $a$ is an element of field and $X$ and $Y$ are $m\times n$ and $n\times p$ matrices
QUESTION [0 upvotes]: Show that $a(XY)=X(aY)$ where $a$ is an element of field and $X$ and $Y$ are $m\times n$ and $n\times p$  matrices.

REPLY [2 votes]: By using the notation to designate the coefficient of the $i$th row and $j$th column of a matrix $A$: $(A)_{i,j}$, we have
$$(a(XY))_{i,j}=a\cdot (XY)_{i,j}.$$
On the other hand,
\begin{align}
(X(aY))_{i,j} & =\sum_{k=1}^{n} (X)_{i,k} (aY)_{k,j}
\\ & =\sum_{k=1}^{n} (X)_{i,k}\  a\cdot(Y)_{k,j}
\\ & = a\cdot\sum_{k=1}^{n} (X)_{i,k} (Y)_{k,j} = a\cdot(XY)_{i,j}.
\end{align}


TITLE: Finding a perpendicular vector from a line to a point
QUESTION [4 upvotes]: Let's say I have an arbitrary linear line described by $y=mx+b$. I also have a point $P(x_1,y_1)$ that is not on that line. I suppose $P$ can be anywhere relative to the line. How can I find the vector $v$ the stretches from the point to the line (and is perpendicular to it)?
Here's what I tried: I figured that the difference in the $x$ and $y$ components ($\delta y$ and $\delta x$) basically describe a vector twice as large as $v$ but in the same direction. So I tried finding $\delta x$ and $\delta y$ by the following:
$\delta y=y_1-(m\cdot x_1 +b)$
$\delta x = x_1-(y_1/m - b/m)$
The vector V is then $[\delta x/2, \delta y/2]$. I divide by two because I figure the vector described by $[\delta x, \delta y]$ goes all the way to the top left corner of the imaginary rectangle that is completed by $\delta x$ and $\delta y$.
Am I correct? If not, why? Regardless if I am, is there an alternative or smarter way to accomplish this? I would appreciate any input here.

REPLY [0 votes]: Let the line $L$ be
$$ax+by+c = 0 \tag{1}$$
Let P be
$$ P = (x_1, y_1) \tag{2}$$
Now, to make calculations easier, shift the coordinate system by the vector given by the position $P$, so that $P$ is the origin of the new coordinate system and the new equation for $L$ is
$$ a(x+x_1) + b(y+y_1) + c = 0 \tag{3}$$
If we let a the constant $s$ be the constant at the end of this equation, $$s = ax_1 + by_1+ c$$ The equation for $(3)$ becomes much shorter,
$$ax + by + s = 0 \tag{4}$$
To calculate the vector that we need, which is from the line to $P$, we need to write the equation for the perpendicular line to the line $L$ and through the point $P$, which in the new coordinate system is given by $$ bx - ay = 0 \tag{5}$$
After a few manipulations we get $ bx = ay $ and $$x = \frac{a}{b} y$$ Now the intersection of $L$ and its perpendicular yields the other point of the vector. Substituting this into $(4)$ yields
$a\frac{a}{b}y + by + s = 0$ then $(a^2 +b^2)y = -bs$, and finally
$$y = \frac{-bs}{a^2 + b^2}$$
And similarly $x$ can be shown to be
$$x = \frac{-as}{a^2 + b^2}$$
But since the vector begins at $(x,y)$ and goes to $P$, it's $$\begin{pmatrix} -x \\ -y \end{pmatrix}$$ So the vector may be written as
$$\begin{pmatrix} a\\ b \end{pmatrix} \times \frac{ax_1 + by_1 + c}{a^2+b^2}$$
If the line is given as $y = mx +b$, we need only to rewrite it as $mx + b - y = 0$ and the new value is
$$\begin{pmatrix} m\\ -1 \end{pmatrix} \times \frac{mx_1 + b -y_1}{m^2+1^2}$$.


TITLE: Representing a for loop with modulus in formal notation
QUESTION [1 upvotes]: I have the following section of code I am writing for research. Basically I need to formally represent the mathematical notation behind a set that follows:
for (int i = 1; i < h; i++)
    if (h % i == 0)
        // i is part of the set

I am not very familiar with formal mathematical notation, and was wondering if anyone can help me express this.
Thanks

REPLY [0 votes]: The set could be specified like this: $S = \{i : (1 \leq i < h) \land (h ~\%~ i = 0)\}$.