\section{Introduction}
Renewable energy sources have received much attention in the recent years as they provide cleaner and cheaper energy as compared to fuel-based sources. The level of the carbon footprint that is produced from energy production plants has reached an alarming level around the globe. Carbon dioxide (CO2) is one of the major gases that emit from production plants, which has contributed significantly to global warming. The high level of CO2 pollution has adversely affected the ozone layer, the environment, and the quality of human lives in general. Moreover, it is estimated that the level of fossil fuel production will decline, and within a few generations from now, the world will exhaust all its reserve of fuel.  Therefore, there is an urgent need to explore renewable energy sources. Wind energy is one of these sources that can be exploited in many areas. It provides cheaper energy production with less required maintenance when compared to the solar-based solutions. Moreover, wind energy can provide power during the day and night, given the availability of the wind, unlike the solar-based energy. These impressive features have made this energy production sector one of the rapidly growing ones. And as the world energy demand is growing each year, it is expected for this sector to grow even more in the coming years.

Variable speed wind turbines rotate with a variable speed that follows the speed of wind. Therefore, they achieve maximum utilization of the aerodynamic power. Most of the used wind turbines around the world are of the variable speed type due to less interruption in production and higher aerodynamic efficiency. However, the random and highly fluctuating wind speed has showed several challenges for their control. Controlling the wind turbine rotor under variable speed has become an active research topic in the recent years, as the wind industry started to adopt more advanced control techniques due to the deficiency of the classical control methods.

Variable speed wind turbine has three regions of operation. In the first region, the wind turbine is just starting up. In the second region, the turbine has a normal operation mode, and it tries to maximize the capture of wind energy. The last region happens with high wind speed, and in this region, the turbine limits its rotor speed to protect its mechanical and electrical components. This work is focused on the second region of operation with the aim to optimize the wind energy extraction in this region through turbine’s rotor speed control.

Many control approaches have been proposed to maximize the energy extraction from the wind turbine. However, the vast majority of these control approaches were linear and used a linearized version of the turbine model \cite{Boukhezzar2007}. These linear methods may succeed in attaining the desired level of power extraction and other set control requirements. Nevertheless, their performance may suffer greatly when subjected to parameter uncertainties, wind vibration, disturbances, and other unexpected external parameters \cite{Grimble1996}. Therefore the research in this area has been shifted toward exploring nonlinear control methods \cite{AlJodah2020a,AlJodah2020b}. A nonlinear control approach was proposed in \cite{Thomsen2007} to decouple disturbance from the speed control of a variable speed wind turbine. Feedback linearization was utilized in this process to transform the nonlinear system into a linear one. The wind fluctuations were reduced using this method, and a better performance was observed compared to LQG control. In \cite{Evangelista2010}, four nonlinear control strategies were tested for variable speed wind turbines. Sliding mode with variation law, super-twisting,  sub-optimal, and twisting were tested in simulation. It was found that all controllers were able to track the high-speed signal with no issues. Moreover, the super twisting was the easiest from the implementation perspective. Adaptive super twisting was developed in \cite{Evangelista2013} to maximize energy extraction. The strategy showed high robustness against uncertainties. An adaptive sliding mode control strategy was proposed in \cite{ Merabet2011}  for the speed tracking problem of variable speed turbines. The sliding surface is used in the update equation of the control gains, which resulted in improving the tracking results. In \cite{Narayana2009}, wind speed forecasting was used in variable speed turbines for maximum power point tracking. The prediction ability of this method improved its response time when compared to the traditional ones.  A fuzzy control system was developed in \cite{Zhang2006} for maximum power extraction. The Fuzzy controller provided fewer power fluctuations when compared with the classical PID. Another fuzzy controller was proposed in \cite{Yao2009} to maxmize power extraction through controlling turbine speed to follow a reference optimum speed. An active disturbance rejection method was proposed in \cite{Zhang2011} to maximize energy extraction via pitch angle control. The effectiveness of the method was proved in simulation with various wind profiles. A torque control method was proposed in \cite{Liao2009}, to control a variable speed wind turbine. The study investigated the utilization of the rotor side and the generator side speed signal for the torque control. It was found that the generator side provides better stability than the other one. A fuzzy control method was proposed in \cite{Boukhezzar2007} to have smooth output power. The method enhanced power regulation when compared with the PID control. Other control methods were proposed to solve control challenges in wind turbines such as second order sliding mode \cite{Zhang2021}, H$_{\infty}$\cite{FarajiNayeh2020},  switching sector sliding mode \cite{Berrada2020}, just to name a few. However, in these methods, the control action is selected to have a value that is higher than the maximum bound of the disturbance signal. This cause high load on the generator, and chattering in the signal for the discontinuous controllers. Therefor, this works reduces the demand on the control action to be higher than the disturbance estimation error instead of the disturbance itself.

In this work, robust control methodology is proposed to make the rotor's speed of a variable speed wind turbine follow a trajectory that maximize power extraction. A PI-based sliding surface is used to have a high level of tracking performance. The stability analysis was conducted and verified using the Lyapunov stability criteria. Moreover, the necessary conditions for stability were found. A simulation study was conducted to verify the proposed strategy, and the results were compared with the classical PID controller. The results confirmed the high level of accuracy and robustness of the proposed approach.

\section{Modeling}
The block diagram of the wind turbine is shown in Fig. \ref{fig:WT_2mass}. The variable speed wind turbine ia made of generator, aeroturbine, and a gearbox. The turbine rotor captures the aerodynamic power provided by \cite{Thomsen2007}
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.8\linewidth]{two_mass_sys.png}
	\caption{Wind turbine two mass model}
	\label{fig:WT_2mass}
	\vspace{-10pt}
\end{figure}
\begin{equation}
J_{aero} = \frac{1}{2} \rho \pi R^2 C_p(\lambda, \beta) v^{3}_{wind}
\end{equation}
where  $R$ is the rotor radius, $\beta$ is the blade pitch angle, $\rho$ is the air density, $v_{wind}$ wind speed, and $\lambda$ is the tip speed ratio and it is given by
\begin{equation}\label{eq:lambda}
\lambda =  \omega_{rot} R/ v_{wind}
\end{equation}
where $\omega_{rot}$ is the rotor's speed. A nonlinear function of $\lambda$ and $\beta$ can describe the power coefficient $C_p$, and it is given by \cite{Slootweg2003}
\begin{equation}
C_p(\lambda, \beta) = \left( \mu_1 \ \varphi(\lambda, \beta) - \mu_2 \beta - \mu_3 \beta^x - \mu_4 \right)e^{-\mu_5 \varphi(\lambda, \beta)}
\end{equation}
\begin{equation}
\varphi(\lambda, \beta) = \frac{\beta^3 + 6 \times 10^{-6}\beta - 3\times 10^{-3} \lambda + 1}{-0.02 \beta^4 + \lambda \beta^3 -0.02 \beta + \lambda}
\end{equation}
where $\mu_1=110.23$, $\mu_2= 0.4234$, $\mu_3= 0.00146$, $\mu_4= 9.636$, $\mu_5= 18.4$, and $x= 2.14$. The maximum power coefficient is achieved when $beta=0$. This can be seen clearly from the 2D shape shown in Fig. \ref{fig:cp_curve}. The aerodynamic power is relative to the rotor torque by
\begin{figure}[h!]
	\centering
	\includegraphics[width=\linewidth]{cp.eps}
	\caption{Wind turbine power coefficient curve.}
	\label{fig:cp_curve}
\end{figure}
\begin{equation}
\tau_{aero} = \frac{P_{aero}} {\omega_{rot}}
\end{equation}
therefore
\begin{equation}
\tau_{aero} = \frac{1}{2 \omega_{rot}} \rho \pi R^2 C_p(\lambda, \beta) v^{3}_{wind}
\end{equation}
The rotor and generator dynamics equations are given by
\begin{equation}\label{eq:Jr}
J_{rot}\dot{\omega}_{rot} = \tau_{aero} - \tau_{low\_speed} - B_{rot}\omega_{rot} - k_{rot} \theta_{rot}
\end{equation}
\begin{equation}
J_{gen}\dot{\omega}_{gen} = \tau_{high\_speed} - \tau_{elec} - B_{gen}\omega_{gen} - k_{gen} \theta_{gen}
\end{equation}
The gear box ratio is given by
\begin{equation}
\eta = \frac{\tau_{low\_speed}}{\tau_{high\_speed}} = \frac{\omega_{gen}}{\omega_{rot}} = \frac{\theta_{gen}}{\theta_{rot}}
\end{equation}
Substituting the gear ratio in the generator equation gives
\begin{multline}\label{eq:Jgen}
\eta^2 J_{gen}\dot{\omega}_{gen} = \tau_{low\_speed} - \eta \tau_{elec} \\ - \eta^2 B_{gen}\omega_{rot} - \eta^2 k_{gen} \theta_{rot}
\end{multline}
Adding Eq. (\ref{eq:Jgen}) and Eq. (\ref{eq:Jr}) results in
\begin{equation}
J_{t}\omega_{rot} = \tau_{aero} - \tau_{gen} - B_{t}\omega_{rot}
\end{equation}
\begin{multline}
(J_{rot} + \eta^2 J_{gen})\dot{\omega}_{rot} = \tau_{aero} - \eta \tau_{elec} \\ - (B_{rot}+ \eta^2 B_{gen}) \omega_{rot} - (k_{rot} - \eta k_{gen} ) \theta_{rot}
\end{multline}
The model can be reduced to a single mass system, therefore the final dynamics equation is given by
\begin{equation}
J_{t}\dot{\omega}_{rot} = \tau_{aero} - \tau_{gen} - B_{t}\omega_{rot} - k_{t}\theta_{rot}
\end{equation}
where $J_{t} = J_{rot} + \eta^2 J_{gen}$, $J_{gen} = \eta \tau_{elec}$, $B_{t} = B_{rot} + \eta^2 B_{gen}$, and $k_{t} = k_{rot} + \eta k_{gen}$. The stiffness has negligible effect on the dynamics. Thus, by ignoring the stiffness and considering the uncertainty in parameters and external disturbances as a lumped disturbance $d$, the final dynamics equation is given by
\begin{equation}
\dot{\omega}_{rot} = \frac{1}{J_{t}}\tau_{aero} - \frac{B_{t}}{J_{t}}\omega_{rot} - \frac{1}{J_{t}} \tau_{gen} - \frac{1}{J_{t}}d
\end{equation}
\section{Control method design}
It has been explained earlier that wind turbines operate in three regions. The current work is concerned with the second region, where the goal is to extract the highest amount of the wind energy. Therefore the control method objective is make the power coefficient $C_p$ at its peak. According to the $C_p$ curves given in Fig. \ref{fig:cp_curve}, it is desired to design a controller that keeps the $C_p$ at its highest possible value, i.e
\begin{equation}
C_p(\lambda^{optimal}, \beta^{optimal}) = C^{optimal}_{p}
\end{equation}
To maximize $C_p$ it is necessary to keep $\lambda$, and $\beta$ at their optimal values. As indicated in Fig. \ref{fig:cp_curve}, the curve is maximized with $\beta^{optimal} = 0$, thus $\beta$ is set to zero in this work. The tip speed is the ratio of rotor to wind speed as was shown in Eq. (\ref{eq:lambda}). To keep it equivalent to its optimal value, it means the ratio should be kept constant. According to Fig. \ref{fig:cp_curve}, the optimal values for $C_p$ and $\lambda$ are given by $\lambda^{optimal}=6.9$, and $ C^{optimal}_{p} = 0.4412$. The wind speed is out of our control, thus that means the fluctuations in the wind speed should be tracked by the rotor to maintain $\lambda$ constant and equivalent to its optimal value. Therefore
\begin{equation}
\omega^{optimal}_{rot} = \omega_{ref} = \lambda^{optimal} v_{wind}/R
\end{equation}

The control methodology will be designed to make the rotor speed follow $\omega^{optimal}_{rot}$, it will be refer to as $\omega_{ref}$ from now on.
Defining the sliding variable as
\begin{equation}
s = k_p\dot{e} + k_i e
\end{equation}
where $e = \omega_{rot} - \omega_{ref}$, differentiating $s$ to get
\begin{equation}
\dot{s} = k_p(\dot{\omega}_{rot}-\dot{\omega}_{ref}) + k_i e
\end{equation}
substituting the rotor dynamics will result in
\begin{multline}
\dot{s} = \frac{k_p}{J_{t}} \tau_{aero} - \frac{k_p B_{t}}{J_{t}}\omega_{rot} - \frac{k_p}{J_{t}}\tau_{gen} - k_p \dot{\omega}_{ref} \\ + \frac{k_p d}{J_{t}} + k_i e
\end{multline}
The stability is investigated by using the following Lyapunov function
\begin{equation}
V_1= \frac{1}{2}s^2
\end{equation}
to have a stable system it is necessary to have a negative definite $\dot{V}_1$ \cite{AlJodah2020c,AlJodah2021}, thus
\begin{multline}
\dot{V}_1 = s \dot{s} = s \left( \frac{k_p}{J_{t}} \tau_{aero} - \frac{k_p B_{t}}{J_{t}}\omega_{rot} - \frac{k_p}{J_{t}}\tau_{gen}\right. \\ \left.- k_p \dot{\omega}_{ref} + \frac{k_p d}{J_{t}} + k_i e\right)
\end{multline}
The control action $\tau_{gen}$ is designed to have two components as follows
\begin{equation}
\tau_{gen} = u_{eq} + u_{sw}
\end{equation}
The equivalent component is designed to maintain $s=0$, thus guarantee the high accuracy of the motion tracking. Therefore it is designed as
\begin{equation}
u_{eq} = \tau_{aero} - B_{t}\omega_{rot} - J_{t} \dot{\omega}_{ref}  + \frac{J_{t}k_i}{k_p}e + \hat{d}
\end{equation}
where $\hat{d}$ is the estimated disturbance. The switching control action is designed to bring the system states to the sliding variable and to reject the system's uncertainties \cite{Ghafarian2020}. It is designed as
\begin{equation}
u_{sw} = \frac{J_{t}}{k_p}(k_1s+k_2sgn(s))
\end{equation}
substituting $\tau_{gen}$ in $\dot{V}_1$ to check the stability condition as follows
\begin{align}
\dot{V}_1 &= s\left(-k_1 s - k_s sgn(s) + \frac{k_p}{J_{t}}\tilde{d}\right)\\
&= -k_1 s^2 - \left(k_2 |s| -   \frac{k_p}{J_{t}}\tilde{d} |s|\right)
\end{align}
therefore
\begin{align}
-k_1 s^2 - \frac{k_p}{J_{t}} |s|\left(\frac{J_{t}k_2}{k_p} - |\tilde{d}|\right) & \le 0 \\
\end{align}
The stability condition is given by
\begin{equation}\label{eq:stab_cond_v1}
\frac{J_{t}k_2}{k_p} > |\tilde{d}|
\end{equation}
This indicates that the system will be stable as long as the disturbance estimation is bounded and the control gain is higher than this bound. In this work, Adaptive Fuzzy Disturbance Observer (AFDO) is designed to estimate the disturbance and hence reduce the demand on the control action to be only higher than the bound of disturbance estimation error rather than the disturbance itself, as evident from the stability condition above.
According to fuzzy basis function, the disturbance can be estimated by
\begin{equation}\label{eq:fuzzy_output_re}
\hat{d}(\bar{\omega}|\hat{\theta})=\hat{\theta}^T\psi(\bar{\omega})
\end{equation}
where $\bar{\omega}=[{\omega}_{1},\ {\omega}_{2}]^T=[{\omega}_{rot},\ \dot{\omega}_{rot}]^T$, $\hat{\theta}$  is an adaptive parameter, $\psi(\bar{\omega})$ is the fuzzy basis function, both of dimension $p$, where $p=m^2$, $m$ is the number of membership functions. The fuzzy basis function can be defined as
\begin{equation}\label{eq:fuzzy_basis}
\psi_{l_1,l_2}(\bar{\omega})=\frac{\prod_{i=1}^{2} \mu_{A_{i}^{l_i}} (\bar{\omega}_{i})} {\sum_{l_1=1}^{m}\sum_{l_2=1}^{m} \left (\prod_{i=1}^{2} \mu_{A_{i}^{l_i}} (\bar{\omega}_{i})\right )}
\end{equation}
where $l_1,\ l_2 = 1 \to m$. The AFDO design is based on the following assumption.

\textbf{\textit{Assumption \cite{AlJodah2020}}}: the optimal  adjustable parameter $\theta^*$ will result in a bounded disturbance estimation that is given by $\displaystyle\sup_{\omega \in M_{\bar{\omega}}}\left | d(\omega)-\hat{d}(\omega|\theta^*)\right|< \varepsilon$, where $M_{\bar{\omega}}$ is a compact set, and $\varepsilon$ is a non-negative scalar.\newline
The AFDO is defined by
\begin{equation}\label{eq:z_dot}
\dot{z}=-\sigma z + p(\bar{\omega},u,\hat{\theta})
\end{equation}
where $z$ is the observer state, $p(\bar{\omega},u,\hat{\theta})=\sigma \omega_2 + \frac{1}{J_t}\tau_{aero} - \frac{B_t}{J_t}\omega_1-\frac{1}{J_t}\tau_{gen}-\frac{1}{J_t}\hat{d}(\bar{\omega}|\hat{\theta})$, and $\sigma$ is a positive constant. The disturbance estimation error is given by
\begin{equation}\label{eq:FDO_dist_observ_err}
\zeta=x_2-z
\end{equation}
Thus, Eq. (\ref{eq:z_dot}) can be rewritten as
\begin{equation}\label{eq:z_dot_2}
\dot{z}=\sigma\zeta + \frac{1}{J_t}\tau_{aero} - \frac{B_t}{J_t}\omega_1-\frac{1}{J_t}\tau_{gen}-\frac{1}{J_t}\hat{d}(\bar{\omega}|\hat{\theta})
\end{equation}
The adjustable $\hat{\theta}$ is adapted to reconstruct $\hat{d}(\bar{\omega}|\hat{\theta})$ which in turn  estimate $d(\bar{\omega})$. Therefore, the error dynamics should be established as follows
\begin{equation}\label{eq:FDO_error_dynamics}
\dot{\zeta}=-\sigma\zeta + \frac{1}{J_t} \left(\hat{d}(\bar{\omega}|\hat{\theta}) - d(\bar{\omega}) \right )
\end{equation}
According to the previously stated assumption and the capacity of fuzzy system to act as a universal
approximator, $d(\bar{\omega})$ can be defined as
\begin{equation}\label{eq:disturbance_fdo}
d(\bar{\omega})= {\theta^*}^T\psi(\bar{\omega}) + \varepsilon
\end{equation}
where $*$ stands for the optimal value, and  $\varepsilon$ is the estimation error. Thus, the following error dynamics can be found
\begin{align}\label{eq:FDO_error_dynamics2}
\dot{\zeta} & = -\sigma\zeta+ \frac{1}{J_t}\hat{\theta}^T\psi(\bar{\omega}) - \frac{1}{J_t}{\theta^*}^T\psi(\bar{\omega})  -\frac{1}{J_t}\varepsilon\nonumber\\
& = -\sigma\zeta+ \frac{1}{J_t}\tilde{{\theta}}^T\psi(\bar{\omega}) + \epsilon
\end{align}
where $\epsilon = -(1/J_t)\varepsilon$, and $\tilde{\theta}^T=\hat{\theta}^T - {\theta^*}^T$. The adaptive law of $\hat{\theta}$ can be found using the following Lyapunov function
\begin{equation}\label{eq:FDO_Lyapunov}
V_2= V_1 + \frac{1}{2}\zeta^2 + \frac{1}{2\gamma}\tilde{\theta}^T\tilde{\theta}
\end{equation}
A stable feedback system can be guaranteed when $\dot{V}_2$ is negative semi-definite. Differentiating $V_2$ provides
\begin{align}\label{eq:FDO_V_dot}
\dot{V}_2 & = \dot{V}_1 + \zeta \dot{\zeta} +\frac{1}{\gamma}\tilde{\theta}^T \dot{\tilde{\theta}}\nonumber \\
& =\dot{V}_1 -\sigma\zeta^2 + \frac{1}{J_t}\zeta\tilde{\theta}^T\psi(\bar{\omega}) + \zeta\epsilon +\frac{1}{\gamma}\tilde{\theta}^T \dot{\tilde{\theta}}\nonumber \\
& = \dot{V}_1 -\sigma\zeta^2 + \tilde{\theta}^T\left( \frac{1}{J_t}\zeta\psi(\bar{\omega})+ \frac{1}{\gamma}\dot{\tilde{\theta}} \right) +  \zeta\epsilon
\end{align}
Defining $\dot{\tilde{\theta}}$ as
\begin{equation}\label{eq:th_tild_adapt_law}
\dot{\tilde{\theta}}= -\frac{\gamma}{J_t}\zeta\psi(\bar{\omega}) = -\bar{\gamma}\zeta\psi(\bar{\omega})
\end{equation}
where $\bar{\gamma}=\gamma/J_t$, and by using this definition, the following adaptive law can be established
\begin{equation}\label{eq:th_adapt_law}
\dot{\hat{\theta}}=\bar{\gamma}\zeta\psi(\bar{\omega})
\end{equation}
Substituting Eq. (\ref{eq:th_tild_adapt_law}) in Eq. (\ref{eq:FDO_V_dot}) will yield
\begin{align}\label{eq:FDO_V_dot_final}
\dot{V}_2 &= \dot{V}_1 -\sigma\zeta^2 +\zeta\epsilon\nonumber\\
& = \dot{V}_1 - \frac{\sigma}{2}\zeta^2 + \frac{1}{2\sigma}\epsilon^2 - \left( \sqrt{\frac{\sigma}{2}}\zeta - \sqrt{\frac{1}{2\sigma}}\epsilon \right)^2\\
& \leq \dot{V}_1  - \frac{\sigma}{2}\zeta^2 + \frac{1}{2\sigma}\epsilon^2\nonumber\\
& \leq - \frac{k_p}{J_{t}} |s|\left(\frac{J_{t}k_2}{k_p} - |\tilde{d}|\right)  - \frac{\sigma}{2} \left( \zeta^2 - \frac{\epsilon^2}{\sigma^2}\right)
\end{align}
For $\dot{V}$ to be negative definite the following condition for $\zeta$ should hold: $ \zeta^2 \geq \epsilon^2/\sigma^2$, along with the condition given in Eq. (\ref{eq:stab_cond_v1}).

It is worth noting here is that the designed approach needs the differentiation of the reference rotor speed. The reference signal is related to the wind speed, and that has high fluctuations. A low pass filter was utilized with the differentiator to reduce the noise amplification during the differentiation process. The switching action is depending on the nonlinear signum function, which, when used, introduces chattering in the control signal. Therefore, it has been replaced with a more smooth function $tanh$.
\begin{figure}[h!]
	\centering
	\includegraphics[width=\linewidth]{sys_bd.pdf}
	\caption{Feedback system block diagram.}
	\label{fig:sys_bd}
	\vspace{-10pt}
\end{figure}
\begin{figure}[h!]
	\centering
	\includegraphics[width=\linewidth]{smc_bd.pdf}
	\caption{Sliding mode controller block diagram.}
	\label{fig:smc_bd}
\end{figure}
\section{Results}
Simulink was used to simulate the wind turbine system and its control approach. The overall feedback system as was designed in Simulink is shown in Fig. \ref{fig:sys_bd}. The control approach is designed as in Fig. \ref{fig:smc_bd}. The signum function was replaced by $tanh(s/0.001)$ to reduce chattering \cite{Ghafarian2020b,AlJodah2018,AlJodah2018b,Das2020,Ghafarian2020c,humaidi2018}.  The $\dot{\omega}_{ref}$ was produced by differentiating the reference speed and filtering the results to reduce the noise. To compare the results with a PID controller a switching block was used to switch between strategies. The model parameters are presented in Table \ref{tab:param}.
\begin{table}[h!]
	\centering
	\caption{Model parameters.}
	\label{tab:param}
	\begin{tabular}{cccc}
		\toprule
		\textbf{Parameter}  & \textbf{Value} & \textbf{Parameter}  & \textbf{Value} \\
		\midrule
		$\rho$ & 1.29 &  $Bt$ & 0\\
		$Jt$ &  1.5 & $C^{optimal}_{p}$ & 0.48\\
		$R$ & 1.26 & $\lambda^{optimal}$ & 6.9\\
		\bottomrule
	\end{tabular}
	\vspace{-10pt}
\end{table}
A wind profile with a speed mean value of 8m/s  was used in the simulation, as shown in Figs. \ref{fig:wind8}. The sliding mode (SM) method was able to track the reference speed with high accuracy when compared with the PID controller as shown in Figs \ref{fig:wr}, and \ref{fig:err}.The mean square error for the PID controller was 15.6675, while for SM was 4.3684, therefore the SM method has improved the tracking results by 72\%. It was noticed that the PID controller has higher overshoot than the SM.

The tip speed tracking of the optimal value is shown in Fig. \ref{fig:lam}. The mean square error of $\lambda$ tracking for the PID controller was 0.4667, while for SM was 0.1194, therefore the SM method has improved the tracking results by 74\%. Again PID showed high overshoot.

The transient response of the control system has seen high spikes for the control action (the generated torque), Fig \ref{fig:T}. This high initial value of the torque is due to the high magnitude of the error, external disturbance, and uncertainty in the parameters. The control action becomes high to overcome all these uncertainties. This high peak of the control signal soon reduces in intensity after a few seconds, and the curve follows the changes in the wind speed. Since the wind speed is not regular and randomly changes its shape, the control signal follows the same behavior. And it is seen as fluctuating or oscillating with the wind changes. The standard deviation of the torque signal for the PID was  7.7306, where for SM was 10.4137. Form these numbers, it can be seen that the SM control action has a higher variation than the PID controller. This is expected as SM needs to change quickly to follow the speed fluctuations. To measure the energy capture efficiency, the following equation is adopted
\begin{equation}
E = \frac{ \int_0^t P_{gen}(\tau) d\tau}{\int_0^t P^{optimal}_{aero}(\tau) d\tau} =  \frac{\int_0^t \tau_{gen}\omega_{gen}}{\int_0^t \frac{1}{2}\rho \pi R^2 C^{optimal}_p v^{3}_{wind}}
\end{equation}
The energy generation is shown in Fig. \ref{fig:P}. The energy capture efficiency for the PID controller was 88.47\%, while for SM was 89.79\%, therefore the SM method has improved the energy production by 1.32\%.
\begin{figure}[!h]
	\centering
	\includegraphics[width=\linewidth]{wind_speed1.eps}
	\caption{8 m/s mean wind speed profile}
	\label{fig:wind8}
	\vspace{-10pt}
\end{figure}
\begin{figure}[!h]
	\centering
	\includegraphics[width=\linewidth]{wr.eps}
	\caption{rotor speed tracking}
	\label{fig:wr}
	\vspace{-10pt}
\end{figure}
\begin{figure}[!h]
	\centering
	\includegraphics[width=\linewidth]{err.eps}
	\caption{rotor speed tracking error}
	\label{fig:err}
	\vspace{-10pt}
\end{figure}
\begin{figure}[!h]
	\centering
	\includegraphics[width=\linewidth]{lam.eps}
	\caption{tip speed $\lambda$}
	\label{fig:lam}
	\vspace{-10pt}
\end{figure}
\begin{figure}[!h]
	\centering
	\includegraphics[width=\linewidth]{T.eps}
	\caption{generator torque}
	\label{fig:T}
	\vspace{-10pt}
\end{figure}
\begin{figure}[!h]
	\centering
	\includegraphics[width=\linewidth]{P.eps}
	\caption{generated power}
	\label{fig:P}
	\vspace{-10pt}
\end{figure}
\section{Conclusion}
In this study, a robust control methodology was designed to optimize the energy extraction in a variable speed wind turbine. The nonlinear model of the wind turbine was developed, and the aerodynamics of the rotor was described as a nonlinear model that depends on the tip and rotor speed. The optimal tip speed and power coefficient were found from the power coefficient curves. These values than later used in the development of the control method to obtain maximum power production. Since the wind speed is out of the design method's control, the rotor speed was used to maintain the tip speed ratio. The sliding mode control method based on an proportional-integral sliding variable and adaptive fuzzy disturbance observer was utilized to maximize the rotor speed tracking accuracy. A simulation study was performed  using Simulink, and the results indicated the high energy generation and tracking ability of the robust control method compared with a PID control. Future studies can be focused on reducing the oscillation in the torque signal.
\bibliographystyle{ieeetr}


\section{Introduction}
Representation learning algorithms typically produce encodings into a Euclidean space. However, if the \textit{manifold hypothesis} (the assumption that the data is well-approximated by a manifold) is taken seriously then a Euclidean target space is unnecessarily big for encoding this manifold: by the Whitney embedding theorem, to embed an $n$-dimensional manifold into a Euclidean space one may need up to $2n$ dimensions for the ambient space. The situation is even worse if the embedding is required to preserve distances (which is obviously important for metric learning), as the Nash embedding theorem shows that in this case the dimension of the ambient space needed is $\mathcal O(n^2)$ \cite{andrews2002notes}. Additionally, such algorithms seldom regularize the encoding space (i.e. encourage it to look like a prior distribution) to ensure that it is well-behaved.

In this work we develop a framework for using a manifold as a target space of deep representation learning algorithms. Following our prior work \citep{korman2018autoencoding}, which uses manifolds as the latent space of an autoencoder, this is done by having the encoding network output chart embeddings and membership probabilities for an atlas of a manifold. In other words, instead of learning a single encoder $f: \mathcal X \to \mathbb R^d$ (where $\mathcal X$ is the input data and $\mathbb R^d$ is the Euclidean encoding space), our technique learns $n$-such encoders together with a scoring function $q: \mathcal X \to [0, 1]^n$ that determines which encoding output to use for a given input. Figure \ref{S1} shows such an example for $\mathcal X$ a circle. For matching the distribution of embeddings to a manifold prior, we introduce a maximal mean discrepancy (MMD) \citep{gretton2012kernel} loss for manifolds.

\begin{figure}[h]
\centering
\begin{subfigure}{.4\textwidth}
\centering
\begin{tikzpicture}

\draw[b, thick] (1, 0) arc (0:120:1);
\draw[r, thick] (-.5, .866025) arc (120:240:1);
\draw[g, thick] (-.5, -.866025) arc (240:360:1);

 \draw[r, ->] (-1.2, 0) -- (-2, 0)
        node[pos=0.5, auto=left] {$\phi_1$};

 \draw[b, ->] (0.6, 1.03923048)  -- (1, 1.73205081)
        node[pos=0.5, auto=left] {$\phi_2$};

 \draw[g, ->] (0.6, -1.03923048)  -- (1, -1.73205081)
        node[pos=0.5, auto=left] {$\phi_3$};

\draw[r, thick] (-2.2, -.5) -- (-2.2, .5);
\draw[b, thick] (0.6669873 , 2.15525589) -- (1.5330127 , 1.65525589);
\draw[g, thick] (0.6669873 , -2.15525589) -- (1.5330127 , -1.65525589);

\end{tikzpicture}
\end{subfigure} \begin{subfigure}{.4\textwidth}
\[
\argmax q(x) = \begin{cases}
1; &x \in \begin{tikzpicture}
\draw[r, thick](-.5, .866025) arc (120:240:.25);
\end{tikzpicture} \\
2; &x \in \begin{tikzpicture}
\draw[b, thick] (1, 0) arc (0:120:.25);
\end{tikzpicture} \\
3; &x \in \begin{tikzpicture}
\draw[g, thick] (-.5, -.866025) arc (240:360:.25);
\end{tikzpicture}
\end{cases}
\]
\end{subfigure}
\caption{Example of a learned atlas (with $d=1$) of a circle. Instead of a single encoder function, we learn $n=3$ many: $\phi_1, \phi_2$ and $\phi_3$, each a map $\mathcal X \to [0, 1]$. For a given point $x\in\mathcal X$, the learned function $q: \mathcal X \to [0,1]^3$ ranks which of the $\phi_i$ to use.}
\label{S1}
\end{figure}

Our framework is flexible enough to apply to a variety of representation and deep metric learning algorithms and in this paper we do experiments with manifold generalizations of SimCLR \cite{chen2020simple} for representation learning and standard triplet training \cite{schroff2015facenet} for metric learning. Our results show that, in low dimensions, replacing the Euclidean target space of these algorithms with a manifold gives a significant improvement in the learned embeddings. Besides being of theoretical interest, low dimensional embeddings have the practical benefits of having faster distance computations, requiring less storage, providing visualizations (in the case of dimensions 2 or 3) and avoiding issues with curse of dimensionality for downstream tasks from the embedding (such as clustering).

In metric learning, we additionally find that the MMD regularization loss alone provides a significant improvement over the baseline in the typical high-dimensional euclidean encoding spaces used for metric learning.

We note that this is extended form of our work \cite{korman2021selfsupervised}, which does not include our metric learning experiments.

Our main contributions are as follows:
\begin{itemize}
\item In section \ref{sec-manifolds} we describe the general procedure for how to adjust typical deep representation and metric learning algorithms to have a manifold as a target space. We define a (semi)-metric on this space and a two-term regularization loss that encourages the encoding distribution to be a uniform and efficient atlas.

\item In section \ref{section-MSimCLR} we show how to adapt SimCLR to encode into a manifold. We validate our approach with experiments on MNIST \cite{lecun2010mnist}, FashionMNIST \cite{xiao2017/online}, and CIFAR10 \cite{krizhevsky2009learning}. Our results show that in small dimensions (we do experiments with dimensions 2, 4, and 8), we get significant improvements over the baseline at the cost of only a few extra \textit{bits} of information (namely the chart number) in the representation space.

\item In section \ref{section-MTriplet} we show how to adapt standard triplet loss training to have a manifold as encoding space. We run experiments on the common benchmark datasets CUB-200 \cite{welinder2010caltech} and Stanford Cars \cite{krause20133d}. As in the case of SimCLR, we also see a significant improvement over the baseline in small dimensions. While, unlike in the SimCLR results, our performance likely does not close the gap enough with the results in high-dimensions to be worth the performance tradeoff in a practical setting, we find these results of theoretical interest and validating the potential of our general framework. Furthermore, we also do an experiment in dimension 64 (a common choice in the literature) and find that keeping the target space Euclidean but supplementing the triplet loss with using the (completely unsupervised) regularizing MMD loss yields a substantial improvement in performance compared to no regularization.
\end{itemize}

\subsection{Related work}
\noindent{\bf Learning an atlas}
The idea of learning a manifold as an atlas was discussed in \cite{pitelis2013learning} where they learn a collection of linear charts via a generalization of principal component analysis. The recent work \cite{kohli2021ldle}, which operates in the non-parametric setting, uses Laplacian eigenmaps to construct local charts and then stitches them together to form a single representation space. Our earlier work \cite{korman2018autoencoding}, which we build upon, learns an atlas as the latent space of an autoencoder. Another difference is that in that work the regularization of the latent space is via adversarial training instead of the MMD loss we use in this paper.

\noindent{\bf Encoding space regularization}
Regularization of a latent space is typically only done in algorithms that have a decoder/generator aspect, such as variational autoencoders \cite{kingma2013auto} and adversarial autoencoders (AAE) \cite{makhzani2015adversarial}. The effectiveness of using an MMD loss in particular on such an encoding space is used in  \cite{tolstikhin2017wasserstein} and inspired our MMD loss function. The work of \cite{grattarola2019adversarial} is of similar spirit to ours, as they extend the latent spaces of adversarial autoencoders to more geometrically interesting ones, namely constant curvature Riemannian manifolds.

\noindent{\bf Metric learning}
A variety of techniques improve upon the foundational triplet training approach to metric learning established in \cite{schroff2015facenet}. These techniques typically adjust the loss function in some way, whether by focusing on triplet sampling methods \cite{schroff2015facenet, harwood2017smart, xuan2020improved}, using proxies as elements of a triplet \cite{movshovitz2017no, teh2020proxynca++} or approximating the triplet loss \cite{do2019theoretically}. There have also been non-triplet based lossed functions \cite{ustinova2016learning, sohn2016improved, liu2017sphereface, zhai2018classification} and ensemble methods \cite{kim2018attention, xuan2018deep}. We refer to \cite{musgrave2020metric} and references therein for a broader overview of current approaches.

Our work, in contrast, comes from the different direction of changing the shape of and regularizing the encoding space. As far as we know there have not been techniques that try to control the embedding distribution (in a similar way to how the latent space of an autoencoder is regularized) in metric learning via an unsupervised regularization term.

Though for simplicity we supplement standard triplet training in this paper, many of the aforementioned algorithms could be supplemented with our technique, and we leave it to future work to investigate the effectiveness of this.

\noindent{\bf Representation learning}
SimCLR, which we experiment with, is just one of many recent advancements in self-supervised representation learning. Other such algorithms include AMDIM \cite{bachman2019learning}, MoCo(v2) \cite{he2019moco, chen2020mocov2}, and BYOL \cite{grill2020bootstrap}. We leave to future work the results of adapting these to have manifold output.

\section{Manifolds as encoding spaces} \label{sec-manifolds}
Loosely, a $d$-manifold is a space $\mathcal M$ that locally looks like $\mathbb R^d$. This notion is made rigorous by requiring the existence of \textit{charts}, which are are homeomorphisms of an open subset of $\mathbb R^d$ to an open subset of $\mathcal M$ such that every point in the manifold is contained in some chart (such a collection of charts is called an \textit{atlas}). Figure \ref{torus-chart} shows an example of a single chart on a 2-manifold.

\input{torus_chart.tex}

\subsection{Representing a dataset as a manifold} \label{data-manifold}
For formally modeling a distribution of data, $\mathcal X$, as a manifold,  we use the same approach as in our earlier work \citep{korman2018autoencoding}. Namely we posit the existence of a latent space $\mathcal Z = [0, 1]^d \times \{1, \ldots, n\}$ of $n$, $d$-dimensional charts with coordinate maps $\psi_i : [0, 1]^d \to \mathcal X$ that forms an atlas of a manifold. We use the uniform distribution as the prior on $\mathcal Z$ and we let $X, Z, J$ denote the random variables on $\mathcal X, [0,1]^d,$ and $\{1, \ldots, n\}$, respectively. We will denote by $\mathbbm 1_y$ the distribution supported at a single point $y$. In our decoder-free setup, we wish to learn:

\begin{enumerate}
\item The inverse mappings of the $\psi_i$, which we denote by $\phi_i : \mathcal X\to [0,1]^d$ and which satisfy
\[
p(z | J=i, X=x) = \mathbbm 1_{\phi_i(x)}.
\]
\item The chart membership function $q = (q_1, \ldots, q_n): \mathcal X \to [0, 1]^n$ defined by
\begin{equation}
q(x) = \left(p(J = 1 \mid X=x), \ldots, p(J = n \mid X=x)\right). \label{q}
\end{equation}
\end{enumerate}
We can then compute the posterior in terms of $q$ and the $\phi_i$ as
\[
p(z, j \mid x) = p(z \mid j, x) p(j \mid x) = q_j(x) \mathbbm 1_{\phi_j(x)}
\]
which gives the prior on $\mathcal Z$ as
\begin{equation}
p(z, j) = \mathbb E_x q_j(x) \mathbbm1_{\phi_j(x)}, \label{pushforward-distribution}
\end{equation}
which we wish to be uniform.

An additional desire is that we have an efficient atlas in the sense that any point $x$ should be in as few charts as possible. Thus while $p(J)$ should be uniform, we want the conditional distributions $p(J \mid X=x)$ to have low entropy: if the distributions $p(J \mid X=x)$  are mostly deterministic then encoding $x$ requires us to only keep the coordinates and chart number for the chart with highest probability for $x$. In other words, for the representation of $x$ at inference time we take
\begin{equation}
x \mapsto (\phi_i(x), i) \in \mathbb R^d \times \{1, \ldots, n\}, ~~ \text{where } i = \argmax_j q_j(x), \label{compressed-rep}
\end{equation}
which has just $\log n$-more bits of information than the Euclidean case.

\subsection{A Maximal Mean Discrepancy Loss for Manifolds}
For a given embedding task, we propose to parameterize the functions $\{\phi_1, \ldots, \phi_n, q\}$ using neural networks and optimize the parameters via gradient descent for a loss function consisting of a task-specific term (e.g. a contrastive loss function) plus a regularization term that encourages the distribution $p(z, j)$ given by (\ref{pushforward-distribution}) to be close to uniform and the discrete distribution $q(x)$ in (\ref{q}) to be close to a deterministic distribution for each $x \in \mathcal X$.

In the Euclidean case (i.e. when $n=1$), there are two popular ways for regularizing the latent space to match a prior distribution: using adversarial training \citep{makhzani2015adversarial, tolstikhin2017wasserstein} or via an MMD \citep{gretton2012kernel} loss \citep{tolstikhin2017wasserstein}. In \citet{korman2018autoencoding} we used an adversarial loss for a manifold latent space but in this work we use an MMD loss due to better training stability and less hyperparameters to tune. If $P$ and $Q$ are two distributions on a common space and $k$ is a reproducing kernel, then the MMD gives a measure of the difference between the distributions, and is defined by
\begin{equation}
\operatorname{MMD}_k(P, Q)^2 = \mathbb E_{y_1, y_2 \sim P\times P} k(y_1, y_2) - 2 \mathbb E_{y_1, y_2 \sim P\times Q} k(y_1, y_2) + \mathbb E_{y_1, y_2 \sim Q \times Q} k(y_1, y_2). \label{mmd-def}
\end{equation}

Let  $p$ and $q$ be given by (\ref{pushforward-distribution}) and (\ref{q}), respectively, $\mathcal U_\mathcal Z$ denote the uniform distribution on $\mathcal Z$, and $\mathcal U_{\mathcal J}$ denote the uniform distribution on $\mathcal J = \{1,\ldots, n\}$. For the loss function encouraging $p$ to be close to $\mathcal U_\mathcal Z$ we will take an approximation of $\operatorname{MMD}_{k_\mathcal Z}(p, \mathcal U_\mathcal Z)^2$ and for the loss function encouraging $q(x)$ to be far from $\mathcal U_{\mathcal J}$, we take $-\mathbb E_x \operatorname{MMD}_{k_\mathcal J}(q(x),  \mathcal U_\mathcal J)$. This will achieve the goal of encouraging $p$ to be uniform and for $q(x)$ to be deterministic.

To get a kernel $k_\mathcal Z$ on $\mathcal Z$ we can start with a kernel  $k_0$ on $[0, 1]^d$ and then define
\[
k_\mathcal Z: \mathcal Z \times \mathcal Z \to \mathbb R, ~ ((z_1, i), (z_2, j)) \mapsto \delta_{ij} k_0(z_1, z_2)
\]
where $\delta_{ii} = 1, \delta_{ij} = 0$ if $i\ne j$.  For $k_0$ we take, as in \citet{tolstikhin2017wasserstein}, the inverse multiquadratics kernel  but pulled back via the sigmoid function\footnote{this avoids having to compute the final sigmoid activation.}: $k_0(x,y)=\frac{d/6}{d/6 + |\sigma^{-1}(x) - \sigma^{-1}(y)|^2}$. We approximate $\operatorname{MMD}_{k_\mathcal Z}(p, \mathcal U_\mathcal Z)^2$ using the U-statistic estimator in \citet{gretton2012kernel}, adjusted to our manifold setting. Explicitly:

\begin{prop} \label{prop-mmd}
Let $p$ be the distribution on $\mathcal Z$ defined by (\ref{pushforward-distribution}) and let $\mathcal U_\mathcal Z$ denote the uniform distribution on $\mathcal Z$. Given a random sample $\{x_1, \ldots,  x_N\}$ of $\mathcal X$ and a random sample $\{w_1, \ldots, w_N\}$ drawn uniformly from $[0,1]^d$, an estimator for $\operatorname{MMD}_{k_\mathcal Z}(p, \mathcal U_\mathcal Z)^2$ is:
\begin{align}
\ell_\mathcal Z (q, \phi_1, \ldots, \phi_n) &= \frac{1}{N(N-1)} \sum_{\substack{j, k = 1 \\ j \ne k}}^N \sum_{i=1}^n q_i(x_j) q_i(x_k) k_0(\phi_i(x_j), \phi_i(x_k)) \label{ellZ-def} \\
& - \frac{2}{n N^2} \sum_{j, k = 1}^N \sum_{i=1}^n q_i(x_j) k_0(\phi_i(x_j), w_k) +  \frac{1}{nN(N-1)} \sum_{\substack{j, k = 1 \\ j \ne k}}^N k_0(w_j, w_k). \nonumber
\end{align}
\end{prop}
See section \ref{prop-proof} for the proof.

For $\mathbb E_x \operatorname{MMD}_{k_\mathcal J}(q(x),  \mathcal U_\mathcal J)$ we take the kernel $k_\mathcal J :  \{1, \ldots, n\} \times \{1, \ldots, n\} \to \mathbb R, (i, j) \mapsto \delta_{ij}$ and define
\begin{align*}
\ell_\mathcal J(q)  &:= -\mathbb E_x \operatorname{MMD}_{k_\mathcal J}(q(x),  \mathcal U_\mathcal J) = -\mathbb E_x \sum_{i=1}^n \left(q_i(x) - \frac{1}{n}\right)^2.
\end{align*}
The total regularization loss is then
\begin{equation}
\ell_{reg}(q, \phi_1, \ldots, \phi_n) = \lambda_1 \ell_\mathcal Z (q, \phi_1, \ldots, \phi_n) + \lambda_2 \ell_\mathcal J(q) \label{reg-loss}
\end{equation}
for some hyperparameters $\lambda_1, \lambda_2 \in [0,\infty)$.

\subsection{Distance function on a manifold}
To define a metric on the manifold encoding space (which is necessary for some representation learning algorithms), we start with a metric $d_0$ on $[0,1]^d$ and extend it to $\mathcal M$ via
\begin{equation}
d_\mathcal M(x, y) = \begin{cases}
\frac{\sum\limits_i q_i(x) q_i(y) d_0(\phi_i(x), \phi_i(y))}{\sum\limits_i q_i(x) q_i(y)}; & \sum\limits_i q_i(x) q_i(y) \ne 0\\
\infty; &\text{else}
\end{cases} \label{manifold-metric}
\end{equation}
We note that technically, $d_M$ is a semi-metric since the triangle inequality may fail.

\subsection{Summary of our framework} \label{sec-framework}
We now describe our general technique of turning a deep representation learning algorithm to one that has a manifold encoding space. The typical setup for such an algorithm, $\texttt A$, is to learn an encoder $f$ mapping the dataset $\mathcal X$ to the space $\mathbb R^d$ via optimizing some loss function $\ell_\texttt A$ defined on embeddings of a mini-batch $\{x_1, \ldots, x_N\} \subset \mathcal X$.

To adjust the algorithm to one that has a manifold as the encoding space, we see from the discussion in \ref{data-manifold} that $f$ should be replaced by a collection of maps $\phi_1,\ldots,\phi_n, q$ where $\phi_i : \mathcal X \to [0, 1]^d$ is the $i$\textsuperscript{th} coordinate map and $q: \mathcal X\to[0,1]^n$ is the chart membership function. In practice, for the functions $\phi_1,\ldots,\phi_n, q$, we take a backbone network and attach $n + 1$ linear heads followed by a sigmoid activation on the first $n$ (which have output in $\mathbb R^d$) and a softmax activation on the head defining $q$ (which has output in $\mathbb R^n$).

In many contrastive representation learning algorithms (such as SimCLR \citep{chen2020simple}, MoCo v2 \citep{chen2020mocov2}, and BYOL \citep{grill2020bootstrap}), the loss $\ell_\texttt A$ is a function of an auxiliary projection head $h: \mathbb R^d \to \mathbb R^{\tilde d}$ applied to the embedding vectors. In these cases, for the manifold version we use $n$-many projection heads (one for each coordinate chart), and then produce a single projection vector for every data point by taking a sum of these individual projection vectors weighted by $q$. We supplement the resulting loss $\ell_\texttt{MA}$ with the regularization loss (\ref{reg-loss}).

At inference and evaluation time, we use the compressed representation (\ref{compressed-rep}). This ensures that the representation is indeed $d$-dimensional and gives a fair evaluation comparison to $\texttt A$, which embeds into a $d$-dimensional Euclidean space. For example, if $q_i(x)$ were the uniform distribution then using the full-representation $(\phi_1(x), \ldots, \phi_n(x))$ instead of the compressed one would essentially be ``cheating" into an $nd$-dimensional Euclidean representation.

We summarize this procedure in Table \ref{summary-tables} and make note of the following special cases:
\begin{enumerate}
\item If $n = 1$ and $\lambda_1 = 0$ then this corresponds to just the algorithm $\texttt A$ except that $f$ is forced to have range $[0, 1]^d$ (e.g. by applying a sigmoid activation at the end).

\item If $n = 1$ and $\lambda_1 > 0$ then the embedding space is still Euclidean like the previous case, but through the regularization loss the distribution of embeddings is encouraged to be close to the uniform distribution on $[0,1]^d$.

\item The case of $d=2$ yields a powerful data visualization by plotting the input data at their embedding coordinates for their most probable chart. We show such visualizations in the appendix \ref{sec-visualizations}.
\end{enumerate}

\begin{table*}[h]
\caption{Summary of our framework for elevating an algorithm $\texttt A$ to have a manifold encoding space.}
\label{summary-tables}
\centering
\begin{small}
\begin{tabular}{c c}
\toprule
{\sc TYPICAL SETUP} & {\sc MANIFOLD VERSION} \\
\midrule

\makecell{learn a neural network encoder \\ $f : \mathcal X \to \mathbb R^d$} & \makecell{learn a neural network encoder with $n + 1$ heads \\ $f = (\phi_1,\ldots,\phi_n, q):\mathcal X\to [0, 1]^d \times \cdots \times [0, 1]^d \times [0, 1]^n$} \\\\

\makecell{euclidean distance function} & \makecell{metric $d_\mathcal M$ as defined in (\ref{manifold-metric})} \\\\

\makecell{auxiliary projection head $h: \mathbb R^d\to \mathbb R^{\tilde d}$ \\ giving projection map $h \circ f: \mathcal X \to \mathbb R^{\tilde d}$.} & \makecell{$n$ auxiliary projection heads $h_1, \ldots, h_n: \mathbb R^d\to \mathbb R^{\tilde d}$ \\ giving projection map $\mathcal X \to \mathbb R^{\tilde d}, ~ x\mapsto \sum_i q_i(x) h_i(\phi_i(x))$.} \\\\

\makecell{on a minibatch $\{x_1,\ldots,x_N\} \subset \mathcal X$ \\ optimize a loss $\ell_{\texttt A}(f(x_1), \ldots, f(x_N))$.} & \makecell{on a minibatch $\{x_1,\ldots,x_N\} \subset \mathcal X$  optimize a loss \\ $\ell_\texttt{MA}(f(x_1), \ldots, f(x_N)) + \ell_{reg}(f(x_1), \ldots, f(x_N))$.} \\\\

\makecell{Representation of $x\in X$ at inference/evaluation: \\ $f(x) \in \mathbb R^d$.} & \makecell{Representation of $x\in X$ at inference/evaluation: \\ $(\phi_i(x), i) \in \mathbb R^d \times \{1, \ldots, n\}$, where $i = \argmax_j q_j(x)$.}

\end{tabular}
\end{small}
\end{table*}

\section{Applications}
\subsection{MSimCLR} \label{section-MSimCLR}
In this section we generalize SimCLR \cite{chen2020simple} to have a manifold as embedding space and denote the resulting technique by \textit{MSimCLR}. We recall that SimCLR trains a neural network encoder $f: \mathcal X \to \mathbb R^d$ using a projection head $h: \mathbb R^d \to \mathbb R^{\tilde d}$. A mini-batch is formed by choosing $N$-images and augmenting in two different ways, producing examples $\{x_1, \ldots, x_{2N}\}$. The loss function over this batch is a function of the projection head output of the embeddings of these images:
\[
\ell_{\text{SimCLR}}(x_1, \ldots, x_n) = c(h(f(x_1)), \ldots, h(f(x_{2n}))),
\]
where $c$ is a contrastive loss. We follow our meta-procedure from the previous section to adjust SimCLR \citep{chen2020simple} to have a manifold as encoding space (with resulting algorithm denoted by \textit{MSimCLR}) and run experiments on MNIST \citep{lecun2010mnist}, FashionMNIST \citep{xiao2017/online}, and CIFAR10 \citep{krizhevsky2009learning} for $n \in \{1, 4, 16, 32\}$ and $d \in \{2, 4, 8\}$.

In \citet{chen2020simple} representations of SimCLR are evaluated based on the accuracy of a linear classifier trained on top of the representation. In our case, since we have a collection of charts, we evaluate our manifold representation by putting a linear classifier on each chart. We report the mean and standard deviation of the accuracy of our models on the hold out test sets in Table \ref{results-table}, from which we see that our method provides a significant performance boost over vanilla SimCLR, especially in dimensions two and four. Section \ref{sec-setup} outlines the details of our experiments, including hyperparameter selection.

\begingroup
\setlength{\tabcolsep}{2.3pt}
\begin{table*}[h]
\caption{Piecewise linear evaluation accuracy (* denotes no convergence)} \label{results-table}
\centering
\begin{scriptsize}
\begin{tabular}{l c c c c c c c c c c}
\toprule
{\sc \tiny METHOD}  & {\sc \tiny \# CHARTS} & \multicolumn{9}{c}{\sc \tiny DATASET} \\
\cmidrule(l){3-11}
& & \multicolumn{3}{c}{\sc \tiny MNIST} & \multicolumn{3}{c}{\sc \tiny Fashion MNIST} & \multicolumn{3}{c}{\sc \tiny CIFAR10} \\
\cmidrule(l){3-5} \cmidrule(l){6-8} \cmidrule(l){9-11}
& & \multicolumn{3}{c}{\sc \tiny ENCODING DIMENSION} & \multicolumn{3}{c}{\sc \tiny ENCODING DIMENSION}  & \multicolumn{3}{c}{\sc \tiny ENCODING DIMENSION} \\
 & & 2 & 4 & 8 & 2 & 4 & 8 & 2 & 4 & 8 \\
\midrule
SimCLR & - & \valstd{15.3}{6.8} & \valstd{75.4}{2.1} & \valstd{94.5}{1.5} & \valstd{39.8}{9.3} & \valstd{62.7}{1.8} & \valstd{79.3}{0.2}  & * & \valstd{66.2}{0.2} & \valstd{79.6}{0.2} \\
MSimCLR & 1 & \valstd{35.4}{12.4} & \valstd{66.2}{0.5} & \valstd{90.5}{2.0} & \valstd{36.1}{10.8} & \valstd{59.1}{4.9} & \valstd{73.5}{3.9} & \valstd{30.2}{10.8} & \valstd{61.1}{2.4} & \valstd{78.0}{1.7} \\
MSimCLR & 4 & \valstd{75.1}{3.8} & \valstd{89.0}{1.7} & \valstd{94.1}{2.6} & \valstd{62.5}{1.2} & \valstd{71.9}{1.5} & \valstd{79.0}{0.3} & \valstd{54.5}{5.8} & \valstd{68.2}{4.5} & \valstd{81.0}{1.9}   \\
MSimCLR & 16 & \valstd{90.4}{3.0} & \valstd{94.3}{3.5} & \valstd{97.1}{0.2} & \valstd{68.5}{2.2} & \valstd{73.9}{1.0} & \valstd{80.8}{0.1} & \valstd{74.2}{0.8} & \valstd{73.0}{2.0} & \valstd{77.3}{0.5}\\
MSimCLR & 32 & \valstd{91.3}{2.8} & \valstd{96.4}{0.2} & \valstd{96.6}{0.0} & \valstd{72.1}{1.0} & \valstd{72.7}{0.5} & \valstd{74.6}{2.2} & \valstd{71.4}{3.5} & \valstd{64.2}{1.1} & \valstd{73.3}{2.7} \\
\bottomrule
\end{tabular}
\end{scriptsize}
\end{table*}
\endgroup

 \pagebreak

\subsection{MTriplet}\label{section-MTriplet}
In this section we apply our framework from section \ref{sec-framework} to adjust standard triplet training \cite{schroff2015facenet} to have a manifold encoding space. For $d_0$ in \ref{manifold-metric} we use $d_0(x, y) = || \sigma^{-1}(x) - \sigma^{-1}(y)||$, where $|| \cdot ||$ is the usual Euclidean norm, which amounts to computing the Euclidean distance before applying the sigmoid activation in the final chart layers. We denote the resulting technique by \textit{MTriplet}.

\subsubsection{Experiments setup}
We run experiments on the CUB-200 and Stanford Cars datasets using PyTorch. In forming the batches we follow \cite{musgrave2020metric} and sample batches that consist of four images across eight different classes. We use the ``batch-all" strategy and the standard margin choice of 0.2. As is common practice, we start with an ImageNet pretrained GoogLeNet architecture. During training we apply a random resized crop to $227\times227$ and a random horizontal flip with probability 0.5. During evaluation we do a resize to $227\times227$. We use the RMSProp optimizer (as suggested in  \cite{musgrave2020metric}) with learning rate $10^{-5}$.

For both datasets the standard practice in the literature is to use the first half of the classes for training and the second half as a test set. For doing hyperparameter selection we split the train half further, taking the first 80\% classes as training and the remaining 20\% for validation. Our grid search covered the following options:  $\lambda_1 \in \{0.1, 1, 5, 10, 20\}$ and $\lambda_1 \in \{0.05, 0.1, 0.2\}$ and we set $d = 2, n = 16$. We optimized for recall at one (R@1) and  chose the hyperparameters that gave the best average rank across the two datasets. This yielded $\lambda_1 = 1$ and $\lambda_2 = 0.1$. Using these parameters we train on the entire train dataset across $d \in \{2, 4, 8\}$ and $n \in \{1, 4, 16, 32\}$ for 300 epochs, repeating each experiment three times, and report the mean recalls at 1, 2, 4, and 8 on the test sets in tables \ref{cub-results} and \ref{cars-results}. For the baseline comparison we did standard, normalized triplet training in $\mathbb R^{d + 1}$ (i.e., for a fair comparison, we make the target space of the baseline one-dimension higher since the embedding vectors get normalized to the unit sphere, which is $d$-dimensional).

\begingroup
\setlength{\tabcolsep}{4.2pt}
\begin{table*}[h!]
\caption{Results on CUB-200} \label{cub-results}
\vskip 0.15in
\centering
\begin{scriptsize}
\begin{tabular}{c c | c c c c | c c c c | c c c c}
\toprule
{\sc method} & {\sc \# charts} & \multicolumn{4}{c|}{\sc encoding dim 2} & \multicolumn{4}{c|}{\sc encoding dim 4} & \multicolumn{4}{c}{\sc encoding dim 8} \\
\cmidrule(l){3-14}
 &  & R@1 & R@2 & R@4 & R@8 & R@1 & R@2 & R@4 & R@8 & R@1 & R@2 & R@4 & R@8 \\
\midrule
Triplet & & $7.65$ & $13.47$ & $22.40$ & $35.07$ & $14.87$ & $23.08$ & $34.72$ & $48.67$ & $26.93$ & $38.04$ & $50.89$ & $64.84$ \\
MTriplet & 1 & $7.03$ & $12.09$ & $20.12$ & $32.35$ & $14.04$ & $22.20$ & $33.67$ & $48.13$ & $25.24$ & $36.76$ & $49.91$ & $64.02$ \\
MTriplet & 4 & $9.96$ & $17.02$ & $27.13$ & $41.10$ & $17.22$ & $26.94$ & $39.03$ & $53.72$ & $28.67$ & $40.37$ & $53.74$ & $67.22$ \\
MTriplet & 16 & $12.13$ & $20.18$ & $31.33$ & $45.13$ & $19.95$ & $30.33$ & $42.98$ & $56.88$ & $29.77$ & $41.44$ & $54.64$ & $67.71$ \\
MTriplet & 32 & $12.92$ & $20.71$ & $31.75$ & $45.02$ & $19.80$ & $30.02$ & $41.74$ & $55.50$ & $29.73$ & $41.24$ & $54.57$ & $67.96$ \\

\bottomrule
\end{tabular}
\end{scriptsize}

\end{table*}
\begin{table*}[h!]
\caption{Results on Stanford Cars} \label{cars-results}
\vskip 0.15in
\centering
\begin{scriptsize}
\begin{tabular}{c c | c c c c | c c c c | c c c c}
\toprule
{\sc method} & {\sc \# charts} & \multicolumn{4}{c|}{\sc encoding dim 2} & \multicolumn{4}{c|}{\sc encoding dim 4} & \multicolumn{4}{c}{\sc encoding dim 8} \\
\cmidrule(l){3-14}
 &  & R@1 & R@2 & R@4 & R@8 & R@1 & R@2 & R@4 & R@8 & R@1 & R@2 & R@4 & R@8 \\
\midrule

Triplet & & $10.92$ & $18.73$ & $30.85$ & $46.81$ & $18.96$ & $29.91$ & $43.81$ & $59.55$ & $33.53$ & $47.20$ & $61.10$ & $74.14$ \\
MTriplet & 1 & $9.62$ & $17.02$ & $28.26$ & $44.29$ & $16.51$ & $27.15$ & $41.20$ & $56.76$ & $32.78$ & $46.24$ & $60.73$ & $74.11$ \\
MTriplet & 4 & $13.14$ & $22.00$ & $35.52$ & $52.34$ & $22.51$ & $34.56$ & $49.49$ & $64.75$ & $38.77$ & $52.54$ & $65.92$ & $77.94$ \\
MTriplet & 16 & $14.07$ & $24.12$ & $38.27$ & $54.61$ & $25.35$ & $37.89$ & $52.46$ & $67.04$ & $38.77$ & $52.43$ & $65.93$ & $77.74$ \\
MTriplet & 32 & $15.87$ & $25.81$ & $39.51$ & $55.50$ & $25.68$ & $38.52$ & $52.80$ & $67.38$ & $38.78$ & $53.05$ & $66.61$ & $78.13$ \\

\bottomrule
\end{tabular}
\end{scriptsize}

\end{table*}
\endgroup

\subsubsection{MMD loss alone helps in high-dimensions} \label{triplet-reg}
As we moved to dimensions higher than eight, we noticed that increasing the number of charts caused diminishing performance but the case of a single chart gave significant improvements over no regularization. In other words, doing standard triplet training for an embedding network $f: \mathcal X \to \mathbb R^d$ with sigmoid activation, metric defined as in the previous section, and adding the regularization loss term (\ref{ellZ-def}), which in this case simplifies to

\begin{align*}
\ell_\mathcal Z(x_1, \ldots, x_n) = \frac{1}{N(N-1)} &\sum_{\substack{j, k = 1 \\ j \ne k}}^N k_0(f(x_j), f(x_k)) - \frac{2}{N^2} \sum_{j, k = 1}^N k_0(f(x_j), w_k) \\
&+  \frac{1}{N(N-1)} \sum_{\substack{j, k = 1 \\ j \ne k}}^N k_0(w_j, w_k),
\end{align*}
where $\{x_1, \ldots, x_n\}$ is a minibatch and $\{w_1, \ldots, w_N\}$ is a random sample drawn uniformly from $[0,1]^d$.

For this experiment we trained 64-dimensional encodings, again using an ImageNet pretrained GoogLeNet backbone. As in the previous section, we weigh the regularization loss by $\lambda_1 = 1$. For baselines we take the standard, normalized triplet training (denoted ``baseline" in the plots) as well as the case of MTriplet with $n=1, \lambda_1 = 0$ (denoted ``unnormalized" in the plots) to control that the performance improvement is due to the regularization and not that we use unnormalized embedding vectors. We train for 1,000 epochs three times in each case and report in figures \ref{cub-reg} and \ref{cars-reg} the evaluation metrics on the test set as a function of epochs trained. For CUB-200, our results are slightly under the state-of-the art results (using 64 dimensional embeddings with a GoogLeNet backbone) reported in \cite{xuan2020improved}, while for Stanford Cars we outperform the other algorithms at R@4 and R@8.

\begingroup
\setlength{\tabcolsep}{-4pt}
\def0{0}
\begin{figure}[h!]
\begin{tabular}{cccc}
\includegraphics[width=.5\textwidth]{plots/CUB_recall_1.pdf} & \includegraphics[width=.5\textwidth]{plots/CUB_recall_2.pdf} \\ \includegraphics[width=.5\textwidth]{plots/CUB_recall_4.pdf} & \includegraphics[width=.5\textwidth]{plots/CUB_recall_8.pdf}
\end{tabular}
\caption{Test set performance vs. number of epochs trained, on the CUB-200 dataset.} \label{cub-reg}
\end{figure}
\begin{figure}[h!]
\begin{tabular}{cc}
\includegraphics[width=.5\textwidth]{plots/cars_recall_1.pdf} & \includegraphics[width=.5\textwidth]{plots/cars_recall_2.pdf} \\\includegraphics[width=.5\textwidth]{plots/cars_recall_4.pdf} & \includegraphics[width=.5\textwidth]{plots/cars_recall_8.pdf}
\end{tabular}
\caption{Test set performance vs. number of epochs trained, on the Stanford Cars dataset.} \label{cars-reg}
\end{figure}

\endgroup

\pagebreak

\section{Conclusion}
We presented a method for adjusting representation learning algorithms to learn an atlas of a manifold instead of mapping into a Euclidean space. This allows, for a given encoding dimension, more interesting geometries of the embedding space. Our experiments with SimCLR and triplet-loss based metric learning training show that for low encoding dimensions our approach gives much more powerful representations over the baseline. Independent of manifold learning, we also showed that triplet loss training is significantly enhanced by adding a regularizing MMD loss function encouraging the embedding space to be uniform. We hope that this work leads to further research into atlas based manifold learning techniques.


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

The cocktail party effect is the phenomenon that the brain focuses one’s auditory attention on specific stimuli while filtering out other stimuli \cite{getzmann2017switching, bronkhorst2000cocktail}. For example, in complex auditory scenes, people can pay attention to the chat partner while ignoring other people or noise in the background \cite{bronkhorst2000cocktail,cherry1953some}. Researchers in different fields, which include psychology, neuroscience, and information science, have made many efforts to analyze and model the brain’s ability of auditory attention \cite{bregman1994auditory,winkler2009modeling,nelken2008processing}. In the artificial intelligence community, with the rapid development of deep learning, more and more methods have been proposed to model the cocktail party effect. According to the formal definition of the problem, these models can be mainly divided into two types, speech separation and speaker extraction \cite{Wang2018c,michelsanti2020overview}. Suppose that there are $N$ sources, $s_1(t)$, $s_2(t)$, …, $s_N(t)$, then the mixture signal $x(t) \in \mathbb{R}^{1 \times T}$ of $N$ sources can be defined as:

\begin{equation}
  \label{eq:mixture}
  x(t)=\sum_{i=1}^{N}{s_i(t)}.
\end{equation}

\begin{figure}[tbp]
  \centering
  \includegraphics[width=\linewidth]{framework.pdf}
  \caption{The framework of WASE. Note that offset cue and voiceprint cue are optional to injected into speaker extraction module.}
  \label{fig:framework}
\end{figure}

The problem of speech separation can be formulated in terms of estimating $N$ sources, while speaker extraction estimating only one source. Since the speech separation aims to extract each source, this leads to two problems, uncertain number of source problem and label permutation problem \cite{yu2017permutation}. Specifically, $N$ may change dynamically, even unknown sometimes, which result in the uncertain number of source problem; the label permutation problem is that how to add supervision signals to correct channels and avoid training failure caused by gradient conflict. Researchers proposed some solutions represented by deep clustering and permutation invariant training \cite{yu2017permutation,hershey2016deep}. In comparison, speaker extraction takes a different path from the prior. Speaker extraction utilizes additional information that can identify the speaker of interest and regards the information as cues to extract the target speaker. Because speaker extraction is only interested in one speaker, it avoids the two problems above. The additional information includes voiceprint (pitch, identity), visual information (lip movement, facial expression), spatial information (azimuth) \cite{michelsanti2020overview,wang2019voicefilter,afouras2019my,gu2019neural}.

Except the above-mentioned cues, the onset signal is also an important cue in the auditory scene analysis \cite{shamma2011temporal,szabo2016computational}. Previous psychological studies have found that infants are more sensitive to the synchronization between the onset times of auditory and visual stimuli, instead of continuous temporal coherence between auditory and visual events \cite{lewkowicz2010infant,desjardins2004integration}. Infants, who have little experience in language, cannot utilize additional temporal and phonetic information provided by visual signal. However, they can take advantage of the onset cue and benefit from when to hear \cite{werner2009effects}.

In this paper, we propose a novel model that can learn When to Attend for Speaker Extraction (dubbed WASE). WASE can detect the onset signal and regard it as important guidance, as shown in Fig.~\ref{fig:framework}. Specifically, we first obtain the speaker voiceprint from the reference voice. And then we use the onset detector to predict the onset of noisy voice intermediate features with the help of speaker voiceprint. Finally onset signal directly filters intermediate features, which inhibits the signal before onset and guides the extraction module to extract the target speaker from a specific onset time.
We propose a multi-task learning strategy to support our model and use a unified training objective that considers both speaker extraction and onset detection tasks. The oracle onset is generated from the voice activity detection (VAD) results of the clean voice. The training objective of onset detection is to minimize the cross-entropy (CE) loss between oracle onset and the predicted onset in intermediate feature space.
We extend onset cue to onset/offset cues, which means not only the start of voice but also the end of voice is predicted. Interestingly, our onset/offset prediction task is actually speaker-dependent voice activity detection (SDVAD) \cite{chen2020end}. Therefore, the onset/offset-based model completes a combination task of speaker extraction and SDVAD.
Cross fusion of multiple cues is an important mechanism of the brain to process related information from complex scenes and form auditory objects. Therefore, we also extend a single cue to multiple cues, i.e., combine onset/offset cues with voiceprint cue.

Our contributions can be concluded as follow:

\begin{enumerate}
\item As far as we know, it is the first time to explicitly model onset cue in speaker extraction problem, which is considered an important cue for the cocktail party effect in both fields of auditory scene analysis and psychology;

\item We extend the onset cue to onset/offset cues. Since our model reconstructs clean voice and predicts the start/end of voice simultaneously, it completes the speaker extraction task and SDVAD task, which have complementary advantages;

\item We further combine the onset/offset cues with the voiceprint cue. The onset/offset cues play a role in the time dimension while the voiceprint plays a role in the feature dimension, respectively, which is beneficial for the integrity of the auditory object.
\end{enumerate}

\section{RELATION TO PRIOR WORK}
\label{sec:related}
Several prior works for speaker extraction have studied various cues about the target speaker, such as voiceprint \cite{wang2019voicefilter,Xu2020SpEx,koizumi2020speech}, lip movement \cite{afouras2019my,Ephrat2018Looking}, facial appearance \cite{chung2020facefilter}, and spatial information \cite{gu2019neural}. Voiceprint models voice characteristics. Facial appearance models the cross-modal relationship between facial appearance and the voice characteristics of the target speaker. Lip movement information has strong coherence with the dynamic changes of voice, which can enhance the voice signal. The spatial information contained in multi-channel acoustic recordings can filter the speaker in a specific location. Among the cues mentioned above, although lip movement contains onset/offset cues, no study has explicitly modeled it. Therefore, our model is the first time to directly focus on the role of onset/offset cues in speaker extraction.

\section{METHOD}
\label{sec:method}
\subsection{Onset Cue}
\label{ssec:onset}
Onset cue is one of the most important cues in the auditory scene analysis. We take advantage of the reference voice to obtain the onset cue from the mixture input. The reference voice can be represented by $r(t) \in \mathbb{R}^{1 \times {T^r}}$, where $T^r$ denotes the length of the reference voice. $r(t)$ is transformed into a $C$-dimensional representation, ${v \in \mathbb{R}^{C \times 1}}$, by voiceprint encoder. The intermediate features of the mixture can be represented by $U \in \mathbb{R}^{C \times L}$, where $L$ denotes the length of intermediate features. The onset detector searches for the onset signal in $U$ with the guidance of $v$. The onset cue can be defined as:
\begin{equation}
  \label{eq:detector}
  o={\mathcal H}(U, v),
\end{equation}
where $\mathcal H(\cdot, \cdot)$ is the onset cue detector, as shown in Fig.~\ref{fig:model}(C-b). $o \in \mathbb{R}^{1 \times L}$ is a 1-dimensional (1-D) vector of the same length with $U$, and the value of each element in vector represents the current position relationship with onset point, i.e., the value larger than threshold denotes after onset while smaller than threshold denotes before onset. The range of value is $[0, 1]$ and the threshold is $0.5$. That is, the single change point across the threshold is the most important information contained in $o$. To offer guidance in extracting the target speaker voice, $o$ is injected to the intermediate features as follow:
\begin{equation}
  \label{eq:onset_inject}
  U'=U \odot o,
\end{equation}
where $\odot$ denotes element-wise multiplication. $o$ is supervised by oracle onset, which is generated from the clean voice. Specifically, the start time of clean voice is detected by VAD, and we use it to generate a binary vector, ${b \in \mathbb{R}^{1 \times T}}$, with a single change point, i.e., $0$ is before the start time and $1$ after it. But there is a problem with the length of the generated binary vector and $o$ is $T$ and $L$, respectively. To make length equal, we downsample the binary vector with the fixed stride, which is equal to the length of the encoder stride. The downsampled results are named oracle onset in this paper. Note that the stride size of the encoder is $1$ ms in our model, which seems short enough and can hardly influence the precision of oracle onset.

\subsection{Extensions of Onset Cue}
\label{ssec:extensions}
\noindent{\bf{Onset/Offset Cues:}}
Intuitively, offset cue can bring additional benefits. Therefore, we extend the onset cue to onset/offset cues. The structure of the onset/offset detector is the same as the onset detector, as shown in Fig.~\ref{fig:model}(C-b). The only thing that changed is the supervision signal. Specifically, we use both the start time and the end time of clean voice to generate $b$, which has two change points. Then we downsample $b$ to obtain oracle onset/offset cues.

\noindent{\bf{Voiceprint Cues:}}
Onset/offset cues filter mixture intermediate features in the time dimension, while voiceprint filters in the feature dimension. Therefore, we manage to take advantage of both cues in one model. Specifically, the voiceprint cue first filter the noisy voice, and then onset/offset cues filter the noisy voice, as shown in Fig.~\ref{fig:model}(C-c). In contrast, we also show the only voiceprint-based modulation method in Fig.~\ref{fig:model}(C-a).
\begin{figure*}[t]
  \centering
  \includegraphics[scale=0.37]{model.pdf}
  \caption{An illustration of our WASE model. (A): The overall architecture of the proposed WASE. (B): Details of TCN block in speaker extraction module. `D-Conv' indicates a depthwise convolution. (C): Three modulated methods based on different cues.}
  \label{fig:model}
\end{figure*}

\subsection{Model Structure}
\label{ssec:model}
The onset cue and its extensions have been described above. Here we introduce the base model in detail. The structure of the base model is adapted based on our previous proposed framework \cite{hao2020unified}, which includes five modules: voiceprint encoder, onset/offset detector, speech encoder, speech decoder, and speaker extraction module, as shown in Fig.~\ref{fig:model}(A). The model has both the reference voice of the target speaker and the noisy voice as inputs. It is optimized by a multi-task learning strategy, which maximizes the training objective of speaker extraction, and minimizes the loss between the predicted onset cue and the oracle onset cue, simultaneously. The main difference between the current model and the previous model lies in two points: 1) we introduce a more complex module to model more cues, instead of a simple voiceprint cue; 2) we replace the previous feature extraction method of reference voice, short-time Fourier transform (STFT), with the time-coding method, i.e., convolution layer.

\noindent{\bf{Voiceprint Encoder:}}
We obtain the reference voice feature by the time-coding method with a 1-D convolutional layer, whose kernel size is $256$ and stride size is $64$, i.e., $32$ ms and $8$ ms in $8$ kHz. The reference voice feature is processed using two-layer bidirectional long short-term memory (LSTM) and one fully-connected layer. Moreover, we use a mean-pooling to squeeze the time dimension of the reference voice feature and finally get the voiceprint vector $v$. Note that the voiceprint encoder weights will be frozen after $150$ epochs to avoid overfitting.

\noindent{\bf{Speech Encoder \& Decoder:}}
The encoder and decoder follow the same structure in \cite{luo2019conv}, which use a 1-D convolutional layer and a 1-D transpose convolutional layer without bias, respectively. The time-coding method has some advantages over the traditional STFT method, including no phase reconstruction problem, trainable weights, and finer coding granularity.

\noindent{\bf{Extraction Module:}}
We use the temporal convolutional network (TCN) block as the main structure in the speaker extraction module \cite{luo2019conv,Lea2016a}. As shown in Fig.~\ref{fig:model}(B), the TCN block contains pointwise convolution ($1\times1$ Conv), parametric rectified linear unit (PReLU), normalization (Norm), and depthwise convolution (D-Conv). The TCN block has two output paths, one of which is fed to the next TCN block and the other is summarized to obtain the final mask. TCN network has three groups, and each group has eight TCN blocks. The dilation factor of D-Conv increases exponentially by $2$ in each group, which make convolutional receptive fields enlarge layer-by-layer and allow the extraction module to capture correlation in a longer time scale.

We multiply specific cues (onset/offset, voiceprint) with the intermediate features of the speaker extraction module, and inject the results back to guide speaker extraction, as shown in Fig.~\ref{fig:model}(C).

\noindent{\bf{Loss Function:}}
We use the multi-task learning strategy to train our model. One task is reconstructing the target speaker voice, which is supervised by the scale-invariant source-to-noise ratio (SI-SNR) between predicted wave $\hat{s}$ and clean wave $s$. The SI-SNR can be defined as:
\begin{align}
  \label{eq:SI-SNR}
  \begin{cases}
  s_{target}=\frac{\left \langle {\hat{s},s} \right \rangle s}{{\left \| s \right \|}^2};\\
  e_{noise}=\hat{s}-s_{target};\\
  SI-SNR=10 \log_{10} {\frac{{\left \| s_{target} \right \|}^2}{{\left \| e_{noise} \right \|}^2}}, \\
  \end{cases}
\end{align}%
where ${\left \langle {s,s} \right \rangle}$ and ${{\left \| s \right \|}^2}$ are the signal power. The other task is the onset/offset detection, which is supervised by CE loss between the predicted onset/offset vector and the oracle onset/offset. The ratio between two losses is 1.

\section{Experiment Settings}
\label{sec:setting}
We performed experiments on speech separation and speaker extraction benchmark dataset WSJ0-2mix \cite{hershey2016deep}. The training set contains $101$ speakers, and the test set includes $18$ unseen speakers. We downsampled all samples to $8$ kHz. For training samples, we selected randomly two speakers and regarded one as the target speaker and the other as the interfering speaker. Then we sampled two voices from the target and regarded them as clean voice and reference voice, respectively. We sampled one voice from the interfering speaker and set it as the interfering voice. The noisy voice is generated by mixing clean voice and interfering voice at random signal-to-noise ratios (SNR) between $-2.5$ dB and $2.5$ dB, and the length is limited to $4$ seconds. Note that we also used interfering voice as a supervision signal to guide interfering voice prediction in the training phase, which is generated by decoding noisy intermediate features masked by residual of the target mask. In the experiments of onset/offset, we pad silence with a random length of $[200, 800]$ ms at the end of the target/interfering voice, which ensure the existence of offset. In the test phase, we modified the test dataset to match the speaker extraction task. Specifically, we set each one of the two speakers in the original test dataset as the target in turn. We obtained $6000$ samples finally, which is twice the number of original test samples.

We use Adam \cite{kingma2014adam} optimizer with an initial learning rate of $1e^{-3}$, decreasing by half if no improvement in recent $10$ epochs. Our model is trained until the performance of the evaluation set is not improved in $10$ consecutive epochs after the learning rate equals $2.5 \cdot e^{-4}$. We built our model in Pytorch.

\section{Results}
\label{sec:results}
To illustrate the effectiveness of the onset/offset cues, we first trained the base model with given oracle onset/offset cues. As shown in Tab.~\ref{tab:oracle cue}, both models achieved good performance, which showed that the onset/offset guidance information is enough to extract the target speaker. It is reasonable that the onset/offset-based model is better than the onset-based model because of more guidance information in onset-offset cues. And then we used the modulated module to predict onset/offset cues and reported metrics of speaker extraction and SDVAD. The signal-to-distortion ratio (SDR) improvements of models with predicted onset/offset cues is not as good as responding oracle onset/offset cues. However, the model based on the predicted onset cue is better than the oracle onset. We analyze this phenomenon and find that the learned onset cue contains some behavior similar to offset cue at the end of target speech, especially when the target speech is much shorter than the interfering speech. In other words, the network spontaneously learns behaviors similar to onset/offset cues. Therefore, onset-only cue leads to better results than oracle onset. We reported the SDVAD results for the frame-level evaluation in terms of accuracy (ACC) and F-score (F1, the harmonic mean of precision and recall). The ACC and F1 results show that our model could capture the target speaker onset/offset signal. The higher layer, the higher ACC and F1 performance, which is more obvious in onset/offset. It indicates that the onset/offset cues of the target speaker are abstract representations, and it is difficult to effective detection only relying on the voiceprint and the representation of low-layer in our model. The visual information may be a better choice for onset/offset detection solely relying on low-layer representation, e.g., lip movement.

\begin{table}[htbp]
  \caption{SDR improvements (dB) with onset/offset cues on WSJ0-2mix dataset.}
  \label{tab:oracle cue}
  \centering
  \begin{tabular}{llll}
    \toprule
    \textbf{Cues used} & \textbf{ACC(\%)} & \textbf{F1(\%)} & \textbf{SDRi(dB)} \\
    \midrule
    oracle onset & - & - & 14.78 \\
    oracle onset / offset & - & - & 17.02 \\
    \midrule
    onset & 91/90/92/87 & 95/94/95/92 & 16.31 \\
    onset / offset & 66/90/93/96 & 78/92/94/97 & 16.75 \\
    \bottomrule
  \end{tabular}
\end{table}

The voiceprint representation is difficult to distinguish people with similar voices, while onset/offset cues are challenging to handle simultaneous voices. The model with all mentioned cues above could combine the advantages of all cues and form a more robust representation of the auditory object. In Tab.~\ref{tab:voiceprint}, compared with the model based on single voiceprint cue, the performance of the model with two cues is improved while the additional parameters are negligible.

\begin{table}[htbp]
  \caption{SDR improvements (dB) with two cues on WSJ0-2mix dataset.}
  \label{tab:voiceprint}
  \centering
  \begin{tabular}{llll}
    \toprule
    \textbf{Cues used} & \textbf{\#Params} & \textbf{SDR(dB)} & \textbf{SDRi(dB)} \\
    \midrule
    voiceprint & 7.5M & 16.27 & 16.14 \\
    onset + voiceprint & 7.5M & 17.12 & 16.99 \\
    onset / offset + voiceprint & 7.5M & 17.18 & 17.05 \\
    \bottomrule
  \end{tabular}
\end{table}

A comparison with the previously published results on the WSJ0-2mix dataset is shown in Tab.~\ref{tab:others}. Although our best performance is slightly worse than SpEx+, our parameter quantity is nearly half of it, which shows the potential of onset/offset cues and multi-cue integration.
\begin{table}[htbp]
  \caption{SDR improvements (dB) with different speaker extraction methods based on WSJ0-2mix dataset.}
  \label{tab:others}
  \centering
  \begin{tabular}{lllll}
    \toprule
    \textbf{Methods} & \textbf{\#Params} & \textbf{SDRi(dB)} \\
    \midrule
    SBF-MTSAL~\cite{Chenglin2019Optimization} & 19.3M & 7.30 \\
    SBF-MTSAL-Concat~\cite{Chenglin2019Optimization} & 8.9M & 8.39 \\
    SpEx~\cite{Xu2020SpEx} & 10.8M & 14.6 \\
    SpEx+~\cite{xu2020spex_plus} & 13.3M & {\bf 17.2} \\
    WASE (onset / offset + voiceprint) & 7.5M & 17.05 \\
    \bottomrule
  \end{tabular}
\end{table}

\section{CONCLUSION}
\label{sec:conclusion}
Inspired by auditory scene analysis and psychology, we focus on the onset cue and verify the effectiveness in speaker extraction. Then we explicitly model onset cue using reference voice and obtain comparable performance on the benchmark dataset. We extend onset cue to onset/offset cues, and the performance further improved. Since the onset/offset-based model completes the combination tasks of speaker extraction and speaker-aware voice activity detection, it verifies the mutual benefits between the two related tasks. We also combine the onset/offset cues and voiceprint cue. Onset/offset cues model start/end time of voice and voiceprint cue models the voice characteristics. Therefore, the combination of these perceptual cues can promote the integrity of auditory objects. The experimental performance is also further improved.

\section{Acknowledgments}
This work was supported by the Major Project for New Generation of AI (Grant No. 2018AAA0100400), and the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB32070000).

\vfill\pagebreak

\bibliographystyle{IEEEbib}