1 Introduction↩︎

The digital transformation, real-time analytics, and Artificial Intelligence (AI) are advancing manufacturing toward interconnected Industry 4.0 ecosystems, but cybersecurity remains deficient, leaving legacy plant floor assets exposed to sophisticated threats [1]–[3]. Notable incidents like WannaCry shutdowns at auto plants and the 2019 LockerGoga attack on Norsk Hydro highlight the vulnerability of Machine Tool Controllers (MTC) in managing Computer Numerical Control (CNC) machinery and other equipment on the plant floor [3]. Compounding this vulnerability, most MTCs have nonlinear, closed architectures, limiting insight into their mechanisms and restricting efforts to understand and mitigate their security risks [4]. Among cyberattacks, replay attacks are especially dangerous as they need no model knowledge; attackers can just record and replay legitimate measurement streams, bypassing intrusion detection and risking part quality and catastrophic damage [5], [6]. A common method to detect replay attacks is using authentication signals like watermarking, which are unknown to attackers [5].

Physical watermarking verifies system integrity and authenticity, similar to how traditional watermarks prevent piracy and confirm ownership. Watermarking embeds unique authentication signals into the system to distinguish legitimate from replayed data [7]. However, the effectiveness of this method depends heavily on the careful design of these signals. High detection accuracy may degrade control performance, as overly sensitive detection mechanisms can disrupt normal controller operation [6], [8]. Static or poorly tuned watermarks hinder performance or fail to address evolving attacks [6], [8], [9]. This tradeoff motivates an adaptive watermarking paradigm, an approach that offers greater flexibility but increases complexity, posing additional challenges for nonlinear systems [10]. A promising approach would use adaptive watermarking to detect replay attacks, balancing detection accuracy and system responsiveness. In this paper, Dynamic Watermarking (DWM) refers to adaptive watermarking (online selection of watermark statistics using system feedback) whereas offline fixed-statistics designs are termed static or stationary watermarking.

1.1 Related Works↩︎

1.1.1 MTC Cybersecurity↩︎

The shift to IoT-enabled, data-driven Industry 4.0 workflows has expanded the cyberattack surface of MTCs [1]–[3]. MTCs are vulnerable due to outdated systems without regular security updates, low operator awareness, and dense network connections [2], [3]. Manufacturing facilities should integrate IoT-specific security, such as multi-layer authentication, tamper-resistant encryption, and real-time surveillance, into their operational technology infrastructure [2]. AI and machine learning algorithms can analyze controller data to detect nuanced anomalies while adjusting defense policies. Recent cybersecurity studies on CNC machines have further explored deep learning-based models for detecting cyber-physical attacks, but these mainly serve as passive anomaly detectors rather than active authentication or watermarking mechanisms [11]. Viewing cybersecurity as a core design element, rather than just an operational expense, is crucial for resilient smart manufacturing and MTC operations [2], [3].

1.1.2 Attacks on MTCs↩︎

Cyberattacks on MTCs can cause not only data breaches but also significant physical damage. Among these, deception attacks are especially dangerous because they exploit the trust between cyber and physical components [5], [10]. In a deception attack, an adversary alters system data to cause harmful actions by the system or its users. Three major types of deception attacks on MTCs [10] commonly discussed in the literature are:

• Flip attacks: Flip attacks compromise data and control signal integrity in industrial control systems and MTCs by reversing their sign. In MTCs, this sign reversal in actuator commands leads to large errors and instability as opposing signals accumulate over time [3].

• False-Data Injection (FDI) attacks: FDI attacks undermine data integrity in Cyber-Physical Systems (CPS) by injecting deceptive data into data streams. This can mislead control mechanisms into harmful actions by exploiting the trust among sensors, controllers, and actuators, creating discrepancies between perceived and actual states [1], [3], [12].

• Replay attacks: In a replay attack, adversaries capture and retransmit valid signals to trick the system without prior system knowledge. They record nominal operational data during normal operation and later replay it. This makes the detector and controller operate on a falsified but statistically consistent measurement stream. Meanwhile, attackers manipulate the actuator to drive the plant’s actual behavior away from correct operation. This masking of the measurement channel is precisely why replay attacks are difficult to detect with conventional integrity- and authenticity-based mechanisms [5].

Replay attacks are particularly challenging because the detector and controller continue to observe a benign-looking measurement stream even while the physical plant evolves under manipulated actuation. By contrast, direct falsification mechanisms such as injection or flip attacks often create immediate inconsistencies in the observed measurements themselves.

1.1.3 Replay Attack Detection↩︎

A common replay-attack detection method uses watermarking, assumed unknown to adversaries [5], [13]. It enables prompt detection by disrupting the watermarking during an attack. Physical watermarking is classified into input-added, which alters input signals, and output-added, which alters output signals [14]. This ensures detection of any interference, regardless of whether inputs or outputs are compromised. Mo et al. [15] introduced static watermarking for replay attack detection in controlled systems by embedding authentication signals into control inputs, enhancing detection while analyzing the resulting detection-performance trade-offs. Later research has divided into three complementary streams:

1. Operator–specified stationary watermarking design: These frameworks show that combining random Gaussian signals can expose tampering even in noisy or partially observed plants, including single- and multi-input-multi-output extensions with robustness to non-Gaussian disturbances [16], [17], Linear Time-Invariant (LTI) to linear time-variant extensions [10], lightweight key-based recursion for unsecured channels [18], and multi-layer industrial deployments that reduce false alarms and localize faults [19]. However, these methods mostly assume LTI–Gaussian models and keep the operator-chosen watermark distribution stationary, leading to a static detectability–performance trade-off.

2. Optimization-based offline or online stationary watermark design: Systematic watermarking strategies balance detection accuracy and control performance using control-theoretic optimization. For LTI systems with full model knowledge, they choose a static watermark covariance offline by minimizing a Linear-Quadratic-Gaussian (LQG) cost or maximizing Kullback–Leibler (KL) divergence, producing static but provably optimal signals for cost or detection metrics [6], [20]–[24]. These approaches, however, assume LTI–Gaussian dynamics and a single offline watermark power, so they cannot adapt to changing plant dynamics and are less effective against replay attacks in time-varying or nonlinear MTCs. For unknown LTI parameters, Liu et al. [8] propose an online identification-and-design loop that updates the watermark covariance during a learning phase and provably converges to the optimal stationary covariance. Yet the steady-state watermark remains stationary and still relies on LTI–Gaussian assumptions, so it cannot provide the dynamic adaptation required for nonlinear or time-varying MTCs.

3. Non-stationary, scheduled, and learning-based watermarking design: Recent works relax the stationary watermarking by changing when or how watermarks are applied: periodic schedules reduce control cost for discontinuous replay attacks [25] and have been extended to application-level periodic watermarking with compensation for unmanned marine vehicles [26]; switching multiplicative schemes watermark sensor outputs through hybrid generator-remover filters with time-varying parameters unknown to the adversary [27]; and event-based physical watermarks use innovation-dependent stochastic or deterministic triggering rules to improve the detection–performance trade-off [9], [28]. These methods show the benefit of non-stationary watermarking, but they are mainly based on predesigned schedules, switching and triggering rules, or known linear CPS models, and generally do not provide a detector-driven adaptation mechanism. More recently, RL has been used to regulate watermark timing and strength for replay detection in linear CPS [29]. However, this approach is developed under the LTI assumption and regulates watermark activation or scalar strength, rather than learning a general watermark-statistics policy for nonlinear or time-varying MTC operation.

1.2 Research Gaps & Contributions↩︎

Despite cybersecurity advancements, detecting replay attacks in MTCs remains challenging. The review of existing literature identifies three unresolved shortcomings that are particularly pronounced in proprietary, closed-architecture MTCs:

1. LTI-Gaussian dependence. Most watermarking frameworks assume stationary LTI dynamics and i.i.d. Gaussian noise. Time-variant dynamics and unmodeled nonlinearities common in modern plant floors and MTCs could undermine detection and performance assurances.

2. Limited adaptive watermark-statistics design. Static covariances cannot adapt to operating changes, while non-stationary designs are mainly based on predesigned schedules, switching and triggering rules, or scalar watermark activation and strength adaptation. This leaves a fragile security-performance trade-off when watermark sensitivity varies across nonlinear or time-varying MTC regimes.

3. Lack of detection uncertainty quantification. Existing watermark-based replay detectors typically make threshold, schedule, or triggering decisions from residual statistics, but do not maintain a principled online measure of detection confidence for closed-loop watermark adaptation. This limits adaptation, especially under time-varying or nonlinear dynamics where the statistic distribution can drift.

We address research gaps G1-G3 by presenting DynaMark, a model-augmented Reinforcement Learning (RL) framework that learns and adapts watermarking statistics and intensity online without requiring an analytic plant transition model in the RL policy optimization. Using readily available measurements (e.g., position) and online detection-confidence, DynaMark iteratively refines a watermarking policy that adjusts DWM covariance to enhance replay-attack detection with least energy and control-performance overhead. The key contributions of this paper are as follows.

1. We formulate the DWM as a Markov Decision Process (MDP) and design an RL-based DWM framework that generates watermarking signals with adaptive intensity for replay attack detection on MTCs. This allows the watermark to adapt to the detector’s confidence about the system state and adjust watermark intensity dynamically to achieve a balance between control performance and detection power.

2. The unique design of the reward function considers control performance, energy consumption, and detection power simultaneously. Together with a bounded action set, this design enables systematic adjustment of the security–performance trade-off while enforcing actuator-safe watermark budgets.

3. To enable DWM under potentially time-varying and nonlinear dynamics, we characterize how replay-induced watermark mismatch alters the detector statistic and develop a Bayesian belief update that quantifies online detection confidence and serves as a compact input to the RL policy. We further establish a closed-loop mean-square boundedness guarantee under mild local conditions and bounded watermark budgets.

4. The effectiveness of the proposed framework is demonstrated on a Digital Emulator (DE) of a closed-architecture MTC, a nonlinear Mass-Spring-Damper (MSD) benchmark, and a real-world physical stepper-motor testbed, where DynaMark achieves faster detection and lower control performance degradation than benchmark watermarking schemes.

The remainder of this paper is organized as follows: Section 2 presents the MTC system dynamics, watermarking scheme, replay attack model, and detection-power analysis. Section 3 details the DynaMark methodology, the RL-based algorithm, and stability results. Section 4 reports experiments on the DE, the MSD, and the physical testbed, including comparisons with static watermarking baselines and existing benchmarks. Finally, it concludes in Section 5.

2 Problem Setup↩︎

In this section, we introduce the MTC system dynamics, the DWM scheme, and the replay attack model. The modeling in this section is used to characterize the detector statistic and belief feedback provided to DynaMark. The RL policy optimization is not supplied with an analytic plant transition model. Fig. 1 shows the flowchart of the machine tool monitoring and control process, and we explain each block in detail.

2.1 System Model↩︎

A sensor network monitors the machine tool, as shown in Fig. 1. Let \(\mathbf{x}_t\in\mathbb{R}^{m}\) be the (latent) physical state of the machine tool at time \(t\), evolving according to the discrete-time stochastic nonlinear dynamics \[\label{eq:state} \mathbf{x}_{t+1} = \mathcal{F}_t(\mathbf{x}_t,\mathbf{u}_t) + \mathbf{w}_t,\tag{1}\] where \(\mathcal{F}_t:\mathbb{R}^{m}\times\mathbb{R}^c\to\mathbb{R}^{m}\) is the state transition dynamics, \(\mathbf{u}_t\in\mathbb{R}^{c}\) is the control input from the control logic, and \(\mathbf{w}_t\) is a zero-mean disturbance capturing unmodeled dynamics and exogenous perturbations. The control logic, actuator, and detector do not directly access \(\mathbf{x}_t\); they only observe the measurement (reported output) \(\mathbf{y}_t\in\mathbb{R}^{n}\) given by \[\label{eq:meas} \mathbf{y}_t = \mathcal{H}(\mathbf{x}_t) + \mathbf{v}_t,\tag{2}\] where \(\mathcal{H}:\mathbb{R}^m\to\mathbb{R}^n\) is the sensing map and \(\mathbf{v}_t\) is measurement noise. The control logic and detector use only \(\mathbf{y}_t\), which forms the cyber interface vulnerable to adversarial manipulation (see § 2.3). We therefore model the system dynamics in the measurement space using \[\label{eq:sys} \mathbf{y}_{t+1} = \mathcal{G}_t(\mathbf{y}_t,\mathbf{u}_t) + \mathbf{e}_t,\tag{3}\] where \(\mathcal{G}_t:\mathbb{R}^{n}\times\mathbb{R}^{c}\to\mathbb{R}^{n}\) is a known (possibly time-varying or nonlinear) transition function, and \(\mathbf{e}_t\in\mathbb{R}^{n}\) aggregates the effects of \(\mathbf{w}_t\), \(\mathbf{v}_t\), and modeling mismatch in measurement coordinates. We assume \(\mathbf{e}_{t}\) is independent zero-mean Gaussian distributed, \(\mathbf{e}_{t}\sim\mathcal{N}(0,Q_t)\), with covariance \(Q_t\in\mathbb{S}_{+}^{n}\) symmetric Positive Semidefinite (PSD). Under honest sensing with direct state readout \(\mathbf{y}_t\equiv \mathbf{x}_t\), we have \(\mathcal{H}(\mathbf{x}_t)=\mathbf{x}_t\) and \(\mathbf{v}_t\equiv 0\), so Eq. 1 implies \(\mathcal{G}_t=\mathcal{F}_t\) and \(\mathbf{w}_t=\mathbf{e}_t\).

A feedback controller computes control actions from observed measurements. We model the control logic as a (possibly nonlinear) function \(f_c:\mathbb{R}^{n}\rightarrow\mathbb{R}^{c}\) parameterized by \(\eta\), i.e., \(\mathbf{u}_t = f_c(\mathbf{y}_{\{t\}};\eta)\), where \(\mathbf{y}_{\{t\}}=[\mathbf{y}_0,\mathbf{y}_1,\dots,\mathbf{y}_t]\) is the output history. We assume the detector forms the one-step-ahead prediction \[\label{eq:pred} \widehat{\mathbf{y}}_{t+1} = \mathcal{G}_t(\mathbf{y}_t,\mathbf{u}_t),\tag{4}\] and the corresponding residuals \[\label{eq:res} \mathbf{r}_{t+1} = \mathbf{y}_{t+1} - \widehat{\mathbf{y}}_{t+1}.\tag{5}\] The detector computes \(\{\mathbf{r}_t\}\) and triggers an alarm using preset confidence levels and statistical metrics (see Fig. 1 and § 2.2). It then decides whether the control logic continues generating signals (normal operation) or stops (suspected attack).

2.2 Dynamic Watermarking Scheme↩︎

The main idea of a physical watermark is to inject a random signal, \(\phi_{t}\), which is called the watermark signal, into the system in Eq. 3 . This signal is used to excite the system and check whether the system responds to the watermark signal in accordance with the dynamics of the system. We consider injecting watermarks into the control actions, i.e., \[\label{eq:wm95scheme} \mathbf{u}^{\phi}_{t} = \mathbf{u}_{t} + \phi_{t}.\tag{6}\] In practice, watermarks are randomly drawn from a predefined distribution, typically Gaussian. Practitioners design this distribution to preserve controller performance and normal operations while improving detection. In DWM, however, the distribution or its parameters are adjusted based on the system state and detection performance. We reserve “DWM” for this adaptive setting. By contrast, offline covariance-design methods that produce a stationary covariance \(U\) are referred to as “static watermarking". In this study, we assume the watermark signals \(\{\phi_{t}\}\) are independent zero-mean Gaussian random variables with covariance \(U_{t}\), which is time-varying and (in the DynaMark setting) chosen adaptively.

With the DWM defined in Eq. 6 , the system dynamics are governed by \(\mathbf{y}_{t} = \mathcal{G}_{t-1}(\mathbf{y}_{t-1},\mathbf{u}^\phi_{t-1}) + \mathbf{e}_{t-1}\), while the corresponding state estimate is obtained as \(\widehat{\mathbf{y}}_{t} = \mathcal{G}_{t-1}(\mathbf{y}_{t-1},\mathbf{u}^\phi_{t-1})\). With residuals defined in Eq. 5 , their sequence inherits the same probabilistic attributes as \(\mathbf{e}_{t-1}\), namely \(\mathbf{r}_{t}\sim \mathcal{N}(0,Q_{t-1})\). Therefore, the probability of obtaining the measurement \(\mathbf{y}_{t}\) is computed as \(f_{\mathbf{y}_{t}}(y) = (2\pi)^{-n/2} \det(Q_{t-1})^{-1/2} \exp\left(-1/2 g_{t}\right),\) where \[\label{eq:test} g_{t} = \mathbf{r}_{t}^\top Q_{t-1}^{-1} \mathbf{r}_{t}.\tag{7}\] It is easy to show that \(g_{t}\sim\chi^{2}_{n}\) with \(n=|\mathbf{y}_{t}|\) degrees of freedom [30]. When this probability is low, it means the system is likely to be subject to a certain anomaly or attack. Therefore, the test for detecting an attack involves checking \(g_t \lessgtr \tilde{g}\) where \(\tilde{g}\) is an appropriate threshold. If \(g_{t}\) exceeds the threshold, the \(\chi^2\) detector will trigger an alarm.

2.3 Residuals Analysis under Replay Attack↩︎

Replay attacks aim to keep the detector statistically consistent with nominal operation by replacing the real-time measurement stream with previously recorded legitimate measurements while manipulating control input. We assume the adversary has access to all sensors and controllers and is capable of \((i)\) intercepting and modifying the control commands \(\mathbf{u}_t^a\) sent to the plant, and \((ii)\) observing and altering all sensor readings (let \(\mathbf{y}_t^{\prime}\) be the modified measurements). The attacker implements a two-stage strategy: \((i)\) record a sufficiently long benign trajectory \(\{(\mathbf{y}_t,\mathbf{u}_t)\}\) under nominal closed-loop operation (i.e., without attacker intervention), then \((ii)\) replay the recorded measurements while applying a chosen actuator input sequence to the plant. Thus, during the attack, the detector and control logic receive the spoofed sequence \(\mathbf{y}_t^{\prime}\) while the plant receives \(\mathbf{u}_t^a\) via the actuators. It is conventional to interpret \(\mathbf{y}'_t\) as the output of the virtual (time-shifted) system dynamics in Eq. 8 . This interpretation does not necessarily imply that the attacker runs a virtual system; rather, it serves as a conceptual construct for analysis. \[\label{eq:replay95dyn} \mathbf{y}'_{t+1}=\mathcal{G}_t(\mathbf{y}'_t,\mathbf{u}'_t)+\mathbf{e}'_t, \;\;\;\;\mathbf{e}'_t\sim\mathcal{N}(0,Q_t),\tag{8}\] We denote the residual under replay by \(\mathbf{r}_{t}^A\), where \(\mathbf{r}_{t}^A:=\mathbf{y}'_t-\widehat{\mathbf{y}}'_t\). To characterize \(\mathbf{r}_{t}^A\), we locally linearize \(\mathcal{G}_t\) around the operating point \((\widehat{\mathbf{y}}_t,\Bar{\mathbf{u}}_t)\), where \(\Bar{\mathbf{u}}_t\) is the control signal when \(\widehat{\mathbf{y}}_t\) is estimated. We assume \(\mathcal{G}_t\) is continuously differentiable for all \(t\) and define the Jacobians \(G_t:=\nabla_y \mathcal{G}_t(\widehat{\mathbf{y}}_t,\Bar{\mathbf{u}}_t)\) and \(H_t:=\nabla_u \mathcal{G}_t(\widehat{\mathbf{y}}_t,\Bar{\mathbf{u}}_t)\). For \((\mathbf{y},\mathbf{u})\) near \((\widehat{\mathbf{y}}_t,\Bar{\mathbf{u}}_t)\), a first-order Taylor expansion gives \(\mathcal{G}_t(\mathbf{y},\mathbf{u})\approx \mathcal{G}_t(\widehat{\mathbf{y}}_t,\Bar{\mathbf{u}}_t)+G_t(\mathbf{y}-\widehat{\mathbf{y}}_t)+H_t(\mathbf{u}-\Bar{\mathbf{u}}_t),\) with higher-order terms negligible over the operating region of interest. Then, locally, \[\label{eq:lin95sys} \mathbf{y}_{t+1} \approx G_t \mathbf{y}_t + H_t \mathbf{u}_t + \tilde{\mathbf{u}}_t + \mathbf{e}_t,\tag{9}\] where \(\tilde{\mathbf{u}}_t = \mathcal{G}_t(\widehat{\mathbf{y}}_t,\Bar{\mathbf{u}}_t) - G_t \widehat{\mathbf{y}}_t - H_t \Bar{\mathbf{u}}_t\).

Without DWM and under replay, \(\mathbf{u}_t\) is the nominal control input computed from \(\mathbf{y}'_t\). Because the attacker wants the fake readings \(\mathbf{y}'_t\) to resemble the normal \(\mathbf{y}_t\), it is natural to pair them with the nominal input sequence \(\{\mathbf{u}'_t\}\) that produced them. Additionally, the adversary chooses a replay segment from the same operating regime as the current trajectory so that, for all \(t\) during replay, \((\mathbf{y}'_t,\mathbf{u}'_t)\) stays near \((\widehat{\mathbf{y}}_t,\bar{\mathbf{u}}_t)\), where the local linearization is accurate and the Jacobians vary negligibly. Thus, under consistent replay, \(\mathbf{u}'_t \approx \mathbf{u}_t\) and the Jacobians along replay are approximately nominal. The detector then computes \(\widehat{\mathbf{y}}'_{t+1}=\mathcal{G}_t(\mathbf{y}'_t,\mathbf{u}_t)\), while the replayed measurement satisfies Eq. 8 . Using Eq. 9 , the replay residual satisfies \(\mathbf{r}_{t+1}^A\approx H_t(\mathbf{u}'_t-\mathbf{u}_t) + \mathbf{e}'_t\approx \mathbf{e}'_t\). Hence, they are statistically indistinguishable from nominal ones, and the replay attack cannot be detected by the test in 7 without DWM.

Under DWM in Eq. 6 , \(\mathbf{y}'_t\) is generated by the nominal controller and includes its own watermark realization \(\phi'_t\sim\mathcal{N}(0,U'_t)\). Since the attacker replays past measurements, they generally cannot match the current watermark realization \(\phi_t\sim\mathcal{N}(0,U_t)\) in the replayed data. The lemma below describes the residual distribution during the replay attack.

The proof of Lemma 1 is deferred to Appendix 5.1.

The proof of Theorem 1 is deferred to Appendix 5.2. Moreover, by Lemma 1, \(\mathbf{m}_{t}=\mathbf{0}\) (hence \(\mathbf{b}_{t}=\mathbf{0}\)), so \(g_{t}^A\) is a central generalized \(\chi^2\) (a weighted sum of independent \(\chi^2\) variables). Theorem 1 shows that DWM improves replay attack detection on MTCs. However, choosing the right watermarking intensity is vital for balancing control performance and detection accuracy. These insights will aid in characterizing the feedback from the detector and creating an intelligent DWM framework to address this trade-off.

3 DynaMark Methodology↩︎

This section outlines the main components of the Markov Decision Process (MDP) model for the DWM problem in MTCs, and the method used to solve it. The DynaMark framework is denoted in Fig. 2. A typical MTC and its sensors form the physical environment. At each sampling interval, the agent observes the system state, including the latest \(\mathbf{y}_t\) and the detector’s belief, and chooses \(U_t\) as an action. The watermark generator draws a signal from the distribution, integrates it with the control input, and transmits it to the machine tool. The detector assesses the gap between predicted and actual outputs, updates its belief of replay attacks, and provides an updated state and a reward (balancing detection confidence, control performance, and energy consumption). An RL algorithm uses entities in the replay buffer to refine the policy, enabling the watermark to dynamically adapt to changing operational contexts. Crucially, this continuous modulation of the covariance renders the watermark non-stationary, effectively preventing adversaries from learning the static statistical properties required to design replay attacks.

3.1 Markov Decision Process Formulation↩︎

Formally, an MDP problem [31], [32] is defined as a tuple \(\mathcal{M}=(\mathcal{S},\mathcal{A},\mathcal{T},\mathcal{R},\rho_{0})\) where \(\mathcal{S}\) denotes the state space, \(\mathcal{A}\) is the action space, \(\mathcal{T}:\mathcal{S} \times \mathcal{A} \to \mathcal{P}(\mathcal{S})\) represents the state transition function that is the probability distribution over states given being in state \(s\in\mathcal{S}\) and taking action \(a\in\mathcal{A}\); \(\mathcal{R}:\mathcal{S}\times\mathcal{A}\times \mathcal{S} \to \mathbb{R}\) is the reward function (or expected immediate reward) that is received after transitioning from state \(s\in\mathcal{S}\) to state \(s^{\prime}\in\mathcal{S}\) by taking action \(a\in\mathcal{A}\); and \(\rho_{0}:\mathcal{S}\to[0,1]\) is the starting state distribution. Next, we detail each component in the DynaMark context.

State Space: Let \(s_{t} \in \mathcal{S}\) be the MDP state at time \(t\). We define \(s_{t}=(\mathbf{y}_{t},d_{t})\) where \(\mathbf{y}_{t}\) is the measurement signal, and \(d_{t}\) denotes the \(\chi^2\) detector’s confidence at time \(t\). The detection confidence, \(d_{t}\), indicates the detector’s belief if the system is under attack, using observed \(\mathbf{y}_{t}\). In particular,\(d_{t}\) indicates the detector’s confidence in an attack when an alarm is triggered. On the other hand, if no alarm sounds, \(d_{t}\) illustrates the certainty that there is no attack. The belief is updated sequentially based on the results of the statistical tests. To formalize \(d_t\), define the random variables as follows: let \(\sigma=1\) represent the occurrence of an actual ongoing replay attack, with a value of 0 otherwise; \(I_{t}=1\) if the detector raises an alarm and 0 otherwise. Subsequently, \(d_t := \Pr(\sigma=1|I_{1:t})\).

Action Space: Using DWM requires careful design of the covariance of the watermarking distribution \(\mathcal{N}(0,U_t)\) since the controller performance degrades due to the need for high accuracy in cyberattack detection. Specifically, if \(U_t\) produces \(\phi_t\) with high intensity, the control performance deteriorates, although this results in a marked improvement in the detection accuracy. Conversely, if \(U_t\) yields \(\phi_t\) with low intensity, detection accuracy decreases, but the controller’s performance remains intact. Another advantage of the RL-based DWM is that, unlike static approaches, the RL agent adaptively modulates \(U_t\) over time, creating a non-stationary profile that significantly increases the complexity of the adversarial learning problem. We define the MDP action space \(\mathcal{A}\) as the set of bounded \(c\times c\) symmetric PSD matrices, specifically, \[\label{eq:action95space} \mathcal{A} = \{U \in \mathbb{S}_+^c : \|U\|_F \le U_{\max} \}.\tag{10}\] Here, \(U_{\max}\) is a design parameter reflecting the physical saturation limits of the MTC actuators and a stability budget. In § 3.2, we show that, under a locally contractive nominal closed loop, bounding \(\mathop{\mathrm{\mathbb{E}}}\left[\|\phi\|^2\right]\) via the defined action set in 10 yields a uniform mean-square bound on the tracking error. Then, observing \(s_t\in\mathcal{S}\), the action at time \(t\), \(a_t\in \mathcal{A}\), is the covariance matrix of the DWM distribution (\(U_t\)).

State Transition Function: At time \(t\), observing \(s_t=(\mathbf{y}_t,d_t)\) and taking action \(a_t:=U_t\), a watermark signal \(\boldsymbol{\phi}_t\sim\mathcal{N}(0,U_t)\) is added to the control input as in Eq. 6 . The plant or simulator then generates \(\mathbf{y}_{t+1}\), and the detector updates \(d_{t+1}\) from the alarm output. The RL agent thus observes sampled transitions \((s_t,a_t,r_t,s_{t+1})\) without access to an analytic MDP transition kernel.

The proof of Theorem 2 appears in Appendix 5.3. Equation ?? defines the recursive detector-belief update for the next state component \(d_t\). Under the nominal model in § 2, this update is implemented by propagating the detector statistic through Eqs. 6 , 3 , 4 , 5 , 7 , and ?? . This model-based characterization enables detector feedback but is not supplied as an analytic transition model for RL optimization.

The proof of Lemma 2 is in Appendix 5.4. Equation ?? gives a model-based Type-II error when the assumptions in § 2 hold. Procedure 3 summarizes this computation: given \(I_t\), it uses the current and replayed watermark covariances to form \(\mathcal{S}_t\), estimate \(\mathcal{H}_{g_t^A}(\widetilde{g})\), compute \(\beta_t\), and update \(d_t\). The RL policy observes only the resulting belief \(d_t\), not the closed-form model used to compute it.

Reward Function: Recall the controller performance degrades due to the need for high accuracy in a replay attack detection that requires a careful choice of the DWM covariance. The goal of adaptively deciding on \(U_t\) is to establish a trade-off between detection accuracy and control performance while reacting to the changes in the environment. Let \(r_t(s,a)\) denote the random reward of being in state \(s\) and taking action \(a\) at time \(t\) that is defined as \[\label{eq:reward} r_t(s,a) = -\omega_{1} \| \phi_{t} \|_{1} - \omega_{2} \| \mathbf{y}_{t+1}^{\star} - \mathbf{y}_{t+1} \|_{2} + \omega_{3} |0.5 - d_{t+1}|.\tag{11}\]

This reward is a heuristic shaping objective to balance energy overhead, control performance, and detection confidence; convergence or optimality are not implied by the reward design alone, and are instead addressed via numerical validation later. The term \(\| \phi_{t} \|_{1}\) is used as an actuation-overhead regularizer that encourages sparse watermark excitations (occasional stronger excitation rather than persistent small excitation). This choice reduces wear or overhead while preserving the ability to increase excitation when detection confidence is low. The term \(\| \mathbf{y}_{t+1}^{\star} - \mathbf{y}_{t+1} \|_{2}\) penalizes deviations from the un-watermarked trajectory, maintaining control performance. This deviation term directly shapes the learning objective to reduce the error metric that appears in the stability bound analyzed in § 3.2. The term \(|0.5 - d_{t+1}|\) increases detection confidence by pushing \(d_{t+1}\) toward either 0 or 1, avoiding uncertain detection outcomes. A confidence value close to 0.5 indicates high uncertainty, which weakens detection effectiveness.

While the reward encourages stable operation by penalizing trajectory deviation and watermark magnitude, stability is formally addressed in § 3.2. There, the stability theorem shows that if the nominal closed loop is locally contractive and the DWM covariance is bounded, then the expected tracking error remains uniformly bounded. Thus, the reward acts as a performance-shaping mechanism, while boundedness follows from the closed-loop property and the action constraint.

3.2 Closed-Loop System Stability under DynaMark↩︎

The watermark generated by DynaMark perturbs the commanded control signal by Eq. 6 , where \(\phi_t\sim\mathcal{N}(0,U_t)\) and \(U_t\in\mathcal{A}\) is selected by the RL policy. This subsection establishes a local mean-square boundedness result for the deviation between the watermarked trajectory and the nominal (un-watermarked) trajectory under nominal model in §2.

Let the tracking error be \(\xi_t := \mathbf{y}_t-\mathbf{y}_t^\star\). The following assumptions capture MTCs’ standard operating conditions around a nominal regime and formalize the notion that watermarking stays within actuator limits.

1. Local smoothness of plant dynamics. The closed-loop dynamics \(\mathcal{G}_t\) is continuously differentiable for all \(t\) in a neighborhood of the nominal operating regime. Consequently, the Jacobian of \(\mathcal{G}_t\) with respect to the control input is uniformly bounded along the nominal trajectory: there exists \(\bar H>0\) such that \(\sup_{t\ge 0}\|H_t^\star\|^2\le \bar H\), where \(H_t^\star:=\nabla_u \mathcal{G}_t(\mathbf{y}_t^\star,\mathbf{u}_t^\star)\). Physically, \(\bar H\) quantifies the worst-case sensitivity of the next measurement to small perturbations in the commanded input within the intended operating region.

2. History-dependent controller with bounded local gains. The MTC control logic \(\mathbf{u}_t=f_c(\mathbf{h}_t;\eta)\) uses a finite measurement history \(\mathbf{h}_t=[\mathbf{y}_t,\mathbf{y}_{t-1},\ldots,\mathbf{y}_{t-w}]^\top\) and is differentiable on a compact neighborhood of the nominal histories. Hence, its history Jacobians are uniformly bounded: there exist constants \(\{\bar K_i\}_{i=0}^w\) such that \(\sup_{t\ge 0}\|K_{t,i}\|^2\le \bar K_i\) for all \(i\), where \(K_{t,i}:=\nabla_{\mathbf{y}_{t-i}} f_c(\mathbf{h}_t^\star;\eta)\). This assumption captures that, in normal operation, the controller does not exhibit unbounded amplification of measurement perturbations across the finite memory window.

3. Finite-variance unmodeled disturbances. The aggregate disturbance term (including differences between nominal and actual noise as well as unmodeled effects) has bounded second moment: \(\mathop{\mathrm{\mathbb{E}}}\|\varrho_t\|^2\le \sigma_\varrho^2\) for all \(t\). This assumes process noise, sensor noise, and modeling mismatch stay within finite bounds during normal operation.

4. Local contractivity of the nominal closed loop. Define \(A_t := G_t^\star + H_t^\star K_{t,0}\), where \(G_t^\star:=\nabla_y\mathcal{G}_t(\mathbf{y}_t^\star,\mathbf{u}_t^\star)\) and \(K_{t,0}\) is the instantaneous history gain. Assume the nominal closed loop is locally contractive in the sense that \(\bar A:=\sup_{t\ge0}\|A_t\|<1\). This condition formalizes that, around the operating trajectory, the baseline feedback loop attenuates small deviations rather than amplifying them.

The proof of Theorem 3 is deferred to Appendix 5.6.

3.3 Policy Optimization Method↩︎

Let \(\pi_\theta:\mathcal{S}\rightarrow\mathcal{A}\) denote a parameterized (deterministic) DWM policy with parameters \(\theta\), where the action \(a_t=\pi_\theta(s_t)\) specifies the covariance \(U_t\in\mathcal{A}\) of the watermarking distribution. For a given policy \(\pi_\theta\), define the discounted return from time \(t\) as \(G_t := \sum_{k=t}^{\infty}\gamma^{k-t} r_k(s,a)\), \(\gamma\in(0,1)\), and, the performance objective as the expected return under the trajectory distribution induced by \(\pi_\theta\) and the environment dynamics, \(J(\theta) := \mathop{\mathrm{\mathbb{E}}}_{\pi_\theta}\left[G_0 \big| a=\pi_\theta(s)\in\mathcal{A},s_0\sim\rho_0\right]\). The optimal policy is any \(\pi_{\theta^\star}\) satisfying \(\theta^\star\in\arg\max_\theta J(\theta)\).

Since the action is continuous and the environment is stochastic due to plant disturbances and measurement noise, we adopt a deterministic actor–critic approach that learns both (i) an actor \(\pi_\theta\) and (ii) two critic networks \(Q_{\psi^1},Q_{\psi^2}:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}\) approximating the action-value function \(Q^{\pi_\theta}(s,a) := \mathop{\mathrm{\mathbb{E}}}_{\pi_\theta}\left[G_t \big| s_t=s,\;a_t=a\right].\) We use a DDPG-based policy optimization method [33] with a clipped double-Q critic target inspired by TD3 [34]. We refer to this implementation as DDPG-CDQ. The method keeps the deterministic structure of DDPG, while using two critics and a conservative minimum target to reduce critic overestimation. The actor is optimized using the deterministic policy gradient \(\nabla_\theta J(\theta) = \mathbb{E}_{s\sim d^{\pi_\theta}} \left[ \nabla_\theta \pi_\theta(s) \nabla_a Q_{\psi^1}(s,a)\big|_{a=\pi_\theta(s)} \right],\) where \(d^{\pi_\theta}\) is the discounted state-visitation distribution induced by \(\pi_\theta\).

The choice of DDPG-CDQ is an implementation choice rather than an algorithmic contribution; other continuous-control RL algorithms, such as PPO, SAC, or full TD3, could be used to solve the same DynaMark MDP problem. We use a DDPG-based deterministic actor because the action is a continuous watermark covariance, and the actor can naturally output a factorized representation of this covariance.

Due to the definition of action space in 10 , we reparameterize \(U_t\) through a learned matrix factor to enforce PSD by construction. Specifically, the actor outputs a lower-triangular matrix \(L_t=L_\theta(s_t)\) and we set \(U_t = \Pi_{\mathcal{A}}\left(L_t L_t^\top\right)\), where \(\Pi_{\mathcal{A}}\) denotes a projection onto \(\mathcal{A}\), ensuring both \(U_t\in\mathbb{S}_+^c\) and \(\|U_t\|_F\le U_{\max}\). This parameterization yields a smooth mapping from policy parameters to feasible watermark covariances while respecting the actuator-safe watermark budget.

Given a transition \((s_t,a_t,r_t,s_{t+1})\) sampled from the replay buffer, the DDPG-CDQ target is \(y_t := r_t + \gamma \min_{i\in\{1,2\}} Q_{\psi^i_{\mathrm{tgt}}} \left(s_{t+1},\pi_{\theta_{\mathrm{tgt}}}(s_{t+1})\right).\) The critic parameters are updated by minimizing the Bellman error \(\mathcal{L}(\psi^i) = \mathbb{E} \left[ \left( Q_{\psi^i}(s_t,a_t)-y_t \right)^2 \right]\), for \(i\in\{1,2\}\). The actor is updated by ascending an estimate of \(J(\theta)\) through the critic gradient above. Target networks \((\theta_{\mathrm{tgt}},\psi^1_{\mathrm{tgt}},\psi^2_{\mathrm{tgt}})\) are updated via Polyak averaging, which mitigates oscillations induced by bootstrapping and function approximation. We emphasize that DDPG-CDQ is not a full TD3 implementation. It adopts the clipped double-Q target, but does not use target-policy smoothing or delayed actor updates. As in DDPG, the exploration is introduced by an Ornstein-Uhlenbeck (OU) process [35] prior to projection onto \(\mathcal{A}\), producing a temporally correlated exploration process while preserving feasibility, i.e., \(U_t = \Pi_{\mathcal{A}}\left( (L_t + \Delta_t)(L_t+\Delta_t)^\top \right)\), where \(\Delta_t\) is a OU process.

4 Evaluation Results↩︎

We evaluate DynaMark across three representative settings to assess detection effectiveness and the security–performance trade-off under linear, nonlinear, and hardware-in-the-loop conditions. Implementation details, training hyperparameters, and code to reproduce the experiments are available at: [https://github.com/navidaftabi/DynaMark]. Additional DDPG training details are provided in Appendix 5.7.

4.1 Baseline Watermarking Policies↩︎

We compare DynaMark against three classes of watermarking baselines to isolate the value of learning an adaptive covariance policy from detector feedback.

• Static-covariance watermarking. This baseline uses a fixed watermark covariance throughout the experiment, independent of the measurement state or detector belief. For case studies 1 and 2, the low- and high-intensity constants are selected from empirical lower and upper quantiles of the DynaMark adaptive covariance; for case 2, we also include the empirical mean covariance as an intermediate static operating point.

• Non-RL belief-adaptive watermarking. This baseline uses the same detector belief as DynaMark but replaces RL with the deterministic rule \(U_t = U_{\min} + (U_{\max}-U_{\min})d_t\). It tests whether directly increasing watermark intensity with attack belief suffices, without optimizing a long-term reward over energy, control performance, and detection confidence.

• LTI-based offline and online watermarking. For the physical stepper-motor case, we also compare against the offline covariance design and online watermark update method in [8]. These baselines rely on an LTI surrogate of the motor dynamics and provide a comparison with model-based watermarking designs under mode-switching, time-varying hardware behavior.

4.2 Evaluation Metrics↩︎

We use the following metrics to evaluate nominal-operation cost and replay-detection effectiveness. When computing a norm over vector quantities with components in different units, each component is first normalized by its reference unit so the resulting metric is dimensionless.

• Watermark energy. The watermarking overhead, \(\|\phi_t\|_1\), measures the instantaneous magnitude of the injected watermark signal; lower values indicate less actuation effort and lower watermarking cost.

• Control-performance degradation. The nominal performance loss, \(\|\mathbf{y}_{t}^{\mathrm{\star}}-\mathbf{y}_{t}\|_2\), measures the instantaneous trajectory deviation; lower values indicate smaller impact on nominal control behavior.

• Cumulative control-performance degradation (CPD). To summarize nominal performance, we report \(\mathrm{CPD}=\frac{1}{N}\sum_{i=1}^{N}\sum_t \|\mathbf{y}_{t}^{\mathrm{\star},(i)}-\mathbf{y}_{t}^{(i)}\|_2\). This metric summarizes the total nominal performance cost over the full evaluation horizon.

• False-alarm control (\(\mathrm{ARL}_0\)). We set \(\tilde{g}=\mathcal{H}_{g_t}^{-1}(1-\alpha)\), giving nominal false-alarm probability \(\alpha\) and the average run-length \(\mathrm{ARL}_0\approx 1/\alpha\). When the empirical statistic does not follow the theoretical distribution, we set \(\tilde{g}\) using the empirical \((1-\alpha)\)-quantile of nominal \(g_t\).

• Detection delay (\(\mathrm{ARL}_1\)). Under attack, detection performance is measured by \(\mathrm{ARL}_1=\frac{1}{N}\sum_{i=1}^{N}(T_{d,i}-\tau)\), where \(T_{d,i}=\min\{t\geq \tau:I_t^{(i)}=1\}\) is the first alarm time in replication \(i\). Smaller \(\mathrm{ARL}_1\) indicates faster detection.

• Cumulative post-onset detection uncertainty (CDU). To summarize post-onset belief evolution, we report \(\mathrm{CDU}=\frac{1}{N}\sum_{i=1}^{N}\sum_{t\geq\tau}(1-d_t^{(i)})\). Lower values indicate that the belief moves toward the attack state more quickly after replay begins.

• Post-onset alarm frequency (AF). We report the fraction of post-onset time steps that trigger alarms, \(\mathrm{AF}=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{T-\tau+1}\sum_{t=\tau}^{T}\mathbb{I}\{I_t^{(i)}=1\}\). Higher values indicate more persistent replay detection after attack onset.

4.3 Case Study 1: DE of the Siemens Sinumerik 828D MTC↩︎

4.3.1 Experiment Setup↩︎

The Siemens Sinumerik 828D MTC (see Fig. 4) connects operators and machinery, handling tool changes and processes [4]. This DE replicates the MTC’s 2-axis motion control with high precision by tuning parameterized transfer functions based on real controller data. In our experiments, the motion was analyzed along the \(y-\)axis. The DE uses a linear dynamic model on the measurement space (see Remark 2) with \(A=1.0\), \(B=0.010\), and \(Q=1.3741\times10^{-13}\). The baseline controller is proportional, \(\mathbf{u}_t = K_p(\bar \mathbf{y}- \mathbf{y}_t)\), with \(K_p=1.0\) and \(\bar \mathbf{y}=0.012\). Each rollout has horizon \(T=1200\). The replay-attack prior is \(\sigma\sim\mathrm{Ber}(0.05)\), with onset \(\tau\mid\sigma=1\sim\mathrm{Geom}(1/T)\), and the attacker uses full-history replay. The reward weights are set as \(\omega_1=\omega_2=0.35\) and \(\omega_3=0.3\). Given the data, Theorem 3 applies to this case study with \(\Bar{A}=0.99\), \(\Bar{H}=10^{-4}\), and \(\Bar{K}=0\). By selecting \(\epsilon=0.01\) and any \(\rho<0.00995\), we obtain \(\nu(\epsilon,\rho)<1\), following the verification procedure detailed in Appendix 5.6.1.

We evaluate DynaMark over 40 replications with \(\alpha=0.005\), corresponding to the target nominal false-alarm level used for this case study. The attack scenario involves the attacker collecting full-system measurements, launching a replay attack \(\tau=600\) by replacing true sensor data with recorded data \(\mathbf{y}_t^{'}\), and manipulating the controller action to the machine as \(\mathbf{u}_t^A = -\mathbf{u}_t\) by the end of the horizon. The control logic \(\mathbf{u}'_t\) uses \(\mathbf{y}'_t\), and the detector monitors \((\mathbf{y}'_t,\mathbf{u}'_t)\), while the machine tool receives \(\mathbf{u}_t^A\). This attack cannot be detected without DWM.

4.3.2 Performance Evaluation↩︎

Figure 5 summarizes DynaMark’s average behavior on the DE under nominal operation and replay attack. Under nominal operation, DynaMark policy remains low, the detector statistic stays below the threshold, and the trajectory closely follows the no-watermark baseline. After attack onset, DynaMark increases the variance, producing a residual shift that drives the detector statistic above the threshold and raises the attack belief to at least \(0.95\) within 15 time steps. These results confirm DynaMark’s capacity to balance detection accuracy and control performance in MTCs.

Figure 6 compares DynaMark with the static-variance and non-RL adaptive baseline watermarking policies defined in § 4.1 under nominal operation and replay attack (\(U_{\min}=U_{Q_{0.01}}=10^{-7}\) and \(U_{\max}= U_{Q_{0.99}}=1.9{\times}10^{-3}\)).

1. Under nominal operation, DynaMark keeps watermark energy and control-performance degradation well below the \(Q_{0.99}\) baseline while maintaining substantially stronger excitation than the \(Q_{0.01}\) case. The non-RL adaptive rule also has low nominal cost because the detector belief remains near zero during normal operation.

2. Under replay, DynaMark achieves \(\mathrm{ARL}_1=1\) in all replications and drives the detector belief to one. The \(Q_{0.99}\) baseline also detects quickly but with higher nominal cost, while the \(Q_{0.01}\) and non-RL adaptive baselines are less reliable or slower because their watermark intensity remains too weak before sufficient belief accumulation.

To map the security–performance frontier, we sweep a range of static-variance values, \(10^{-8},10^{-7},\dots,10^{-1},1.0\), and compute \(\mathrm{CDU}\) under attack and \(\mathrm{CPD}\) under nominal operation. Figure 6 (d) shows the expected knee: \(\mathrm{CDU}\) decreases rapidly as watermark intensity increases, with limited performance loss up to the knee, after which additional detection gains level off while \(\mathrm{CPD}\) increases sharply. DynaMark aligns near this knee, indicating that it balances detection and performance without relying on a preset static-variance.

4.3.3 Ablation and Sensitivity Analysis↩︎

We conduct ablation and sensitivity studies on the DE case to assess the contribution of each reward component and the robustness of DynaMark to key design parameters. For the reward ablation, we train separate policies by selectively enabling or disabling the energy, performance, and belief terms in the reward. For the watermark-budget sensitivity, we train separate policies with different values of \(U_{\max}\) while keeping the full reward weights fixed. In contrast, the attack-prior and attack-onset-prior studies are evaluation-only: we use the trained DynaMark policy and vary only \(q_{\mathrm{prior}}\) or \(p_{\mathrm{geom}}\) in the online belief update, while keeping the actual replay attack onset fixed at \(\tau=600\).

Fig. 7 visualizes the aggregated trade-offs between nominal cost and detection effectiveness. The DynaMark policy provides a balanced operating point: it achieves substantially stronger detection than energy-only and no-belief variants, while avoiding the large energy and performance costs of belief-only and no-energy variants. Table ¿tbl:tab:ablation? reports the corresponding aggregated metrics, averaged over evaluation replications and summarized as mean \(\pm\) standard error across seeds. The results show that increasing \(U_{\max}\) improves detection persistence but increases nominal cost, while the evaluation-only prior studies indicate that DynaMark remains effective under moderate prior misspecification, with more conservative priors producing slower belief growth and larger onset-prior values reducing post-onset detection confidence.

4.4 Case Study 2: Nonlinear Mass-Spring-Damper (MSD)↩︎

4.4.1 Experiment Setup↩︎

To assess DynaMark beyond linear dynamics and Gaussian uncertainty, we consider a single-axis nonlinear mass-spring-damper (MSD) system, a canonical reduced-order surrogate for compliant MTC axes and transmission dynamics where structural flexibility and damping shape tracking performance and vibration responses, e.g., axis drive compliance, guideway friction, and coupling elasticity [36]. Such second-order models are widely used in MTC and monitoring because they capture the dominant low-frequency mechanics that interact with feedback loops and are directly exposed to cyber–physical manipulation of sensor feedback and commanded inputs. We adopt the nonlinear dynamics \(m\ddot p + b(\dot{p}) + k(p) = F(t) + u(t)\) with damping (\(b(v)=b_1v+b_2v^3\)) and spring (\(k(p)=k_1p+k_2p^3\)) forces, where \(m\) is the mass, \(p\) is displacement from equilibrium position, \(v=\dot{p}\), \(F(t)\) is a measurable disturbance, and \(u(t)\) is the commanded input (see [37] for more details). With sampling time \(T_s=\Delta t\) and honest sensing, we instantiate Eq. 3 with \(\mathbf{y}_t=[p_t,v_t]^\top\), \(\mathbf{u}_t=[0, u_t+F_t]^\top\), and \(\mathcal{G}_t(\mathbf{y}_t,\mathbf{u}_t)=\left[p_t + T_s v_t;\;v_t - \frac{T_s}{m}\left(b(v_t) + k(p_t) -u_t-F_t\right)\right]\). To test robustness beyond the Gaussian setting used in Lemma 2, we model the uncertainty term \(\mathbf{e}_t\) using a zero-mean heavy-tailed Student-\(t\) distribution with 5 degrees of freedom scaled to match covariance \(Q\). We set \(m=1.0\) kg, \(k_1=0.5\;N/m\), \(k_2=1.0\; N/m^3\), \(b_1=1.0\; N/(m/s)\), \(b_2=0.1\; N/(m/s)^3\), \(T_s=0.01\; s\), and \(Q=10^{-7}I\). The exogenous disturbance is constant (\(F_t=F_0=2.0 \; N\)). Each rollout has horizon \(T=4000\). The replay attack’s prior and onset distribution, and its mechanism, follow the MTC threat model in § 4.3 with \(\tau=1000\). The nominal command is a chirp excitation to induce time-varying operating conditions, i.e., \(u_t^\star=0.1\sin\left(t \omega_t\right) \; N\), \(\omega_t=0.1(\frac{1+F_0}{0.1})^{\frac{t}{T\times T_s}}\). The reward weights are set as \(\omega_1=0.1,\omega_2=0.35\), and \(\omega_3=0.55\).

This case study is intended as an empirical stress test beyond the sufficient closed-loop condition in Theorem 3 and beyond the Gaussian residual setting used in the analytical detector belief characterization. Since the nominal command is programmed rather than measurement-feedback, we do not invoke Theorem 3 here. Instead, stability is enforced through the bounded watermark action constraint and verified by bounded trajectories over all training and evaluation rollouts. Because the residuals are non-Gaussian in this experiment, the miss-detection probability \(\beta_t\) is not computed using the Gaussian closed form in Lemma 2; instead, it is estimated offline by Monte Carlo calibration over finite grids of candidate watermark variances and replay-onset times, and then used during training via a lookup table; see Appendix 5.5 for details.

4.4.2 Performance Evaluation↩︎

Figure 8 shows DynaMark’s behavior on the nonlinear MSD under nominal operation and a replay attack starting at \(\tau=1000\). In nominal operation, the learned policy keeps the DWM variance low, the trajectory remains close to the nominal response, and the detector statistic stays below the calibrated alarm threshold, keeping the posterior attack belief near zero. After replay begins, DynaMark increases the variance, producing a clear residual shift; consequently, the test statistic rises above the threshold and the attack belief quickly saturates to one.

To contextualize the trade-off, Fig. 9 benchmarks against static-variance and non-RL adaptive baselines (\(U_{\min}=U_{Q_{0.05}}=5.7e-7\), \(U_{\text{mean}}=9e-4\), \(U_{\max}=U_{Q_{0.95}}=7e-3\)). As expected, the \(Q_{0.95}\) baseline improves detection performance by strengthening excitation, but it incurs higher watermark energy and control-performance degradation; conversely, the \(Q_{0.05}\) baseline has low nominal cost but delayed and unreliable detection. The mean baseline provides an intermediate static trade-off, reducing \(\mathrm{ARL}_1\) relative to \(Q_{0.05}\) but with higher nominal cost than DynaMark. The non-RL adaptive baseline keeps nominal cost low because the belief remains near zero in normal operation, but under attack its detection delay is more variable, indicating that directly scaling \(U_t\) with \(d_t\) is less reliable than learning a systematic DWM policy.

4.5 Case Study 3: Physical Testbed Implementation↩︎

This section evaluates DynaMark on a physical smart stepper-motor testbed representative of MTC-controlled motion systems used in drilling, milling, laser cutting, and 3D printing. Although the motor plant is simple, the implementation preserves practical controller artifacts such as discrete command updates, communication delays, buffering, and timing jitter, all of which affect replay detection and watermark injection. The firmware and motion-control stack is also closely related to those used in 3D-printer controllers, supporting the relevance of this testbed to manufacturing control pipelines.

4.5.1 Experiment Setup↩︎

The experimental platform integrates a \(\rm NEMA17\) closed-loop stepper-motor with magnetic encoder, driven by an \(\rm MKS\) \(\rm Gen\) \(\rm L\) \(\rm V2.1\) control board and an \(\rm MKS\) \(\rm Servo\) \(\rm 42A\) smart driver as illustrated in Fig 10. The \(\rm MKS\) \(\rm Servo\) \(\rm 42A\)’s firmware allows for servo-like control logic, which was then modified to accept \(U_{t}\) updates via serial commands during motion using a command “watermark \(U_{t}\)”, while ensuring synchronization with the encoder’s closed-loop feedback. This is accomplished via serial communication with a PC running \(\rm MATLAB\) at \(115200\) baudrate with \(\rm CR+LF\) termination. Additional firmware modifications enabled scripted cyberattack scenarios triggered via their respective commands. At each control tick, the driver emits a timestamped state vector \(\langle \boldsymbol{\tau}_t, \mathbf{y}_t, \mathbf{u}_t, \boldsymbol{\phi}_t, \mathrm{attack}_{flag} \rangle\), where \(\boldsymbol{\tau}_t\) are time steps, \(\mathbf{y}_t\) are position measurements (encoder counts converted to mm), \(\mathbf{u}_t\) is the control signal, \(\boldsymbol{\phi}_t\) is the watermark signal, and \(\mathrm{attack}_{flag}\) is a 1/0 variable indicating an active attack.

DynaMark deployment on the motor involves three steps: system identification, DE construction for RL training, and attack monitoring on the physical testbed. For identification, we model the motor as a piecewise linear system over operating points defined by speed and direction. At each operating point, the parameters \(A\), \(B\), \(Q\), and desired set-point \(\overline{\mathbf{y}}\) are estimated by fitting an ARX(1,1) model to watermark-free data \(\{\mathbf{y}_t,\mathbf{u}_t\}\) collected at the firmware’s native acquisition rate of \(\rm 1 \;kHz\). The resulting local models are combined to represent the motor’s nonlinear, mode-dependent behavior. Table 1 reports the identified parameters for the four command segments, and Fig. 11 (a) shows the corresponding continuous motion profile with segment boundaries.

For the physical testbed, Theorem 3 can be verified using the identified local ARX(1,1) models and the implemented firmware controller following steps in Appendix 5.6.1. Although the firmware has a PID structure, the runtime integer-scaled gains are \(K_p=921\), \(K_i=0\), and \(K_d=10\) with scaling factor \(1024\); hence, the implemented controller is a finite-history PD controller with \(w=1\) and local measurement gains \(K_{t,0}=-(921+10)/1024=-0.9092\) and \(K_{t,1}=10/1024=0.0098\). Using Table 1, we obtain \(\Bar{H}=\max_j B_j^2=1.1664\times 10^{-4},\) \(\Bar{K}=K_{t,1}^2=9.54\times 10^{-5},\) and \(\Bar{A}=\max_j |A_j+B_jK_{t,0}|=0.9932\). These values satisfy conditions \((C_1)\) and \((C_2)\). Choosing \(\epsilon=0.005\) gives \(\epsilon<\Bar{A}^{-2}-1\), \(c_\epsilon=6.71\times10^{-6}\), and \(\rho_\epsilon^\star=0.00259\). With \(\rho=\rho_\epsilon^\star\), we obtain \(g_\epsilon(\rho_\epsilon^\star)<\delta_\epsilon\) and \(\nu(\epsilon,\rho)=0.9965<1\).

During online decision-making or attack monitoring, trajectory updates are issued in discrete batches rather than as a fully continuous stream. This discretization leads to stepwise profiles in the recorded motion, even when the commanded trajectory is smooth. Figure 11 (b) illustrates the stepwise behavior from batchwise collection and processing. Nevertheless, the stepper-motor maintains accurate tracking across the entire motion profile under these normal operating conditions. To mitigate the RL agent training cost on the physical testbed, we construct the DE of the motor using these same stepwise trajectories and trained the agent on this DE with the reward weights as \(\omega_1=\omega_2=0.35\) and \(\omega_3=0.3\). This ensures the DE accurately reflects the operational characteristics observed during online decision-making. The last task is performing attack monitoring on the actual testbed. A replay attack scenario is considered wherein the adversary replays the measurements while flipping the control signal. For this physical implementation, the nominal alarm threshold is empirically calibrated from residual data with \(\alpha=0.005\).

The online computational cost is dominated by the belief update. For residual dimension \(d\) and \(K\) Monte-Carlo samples used in replay-side \(\beta_t\) estimation, each detector update costs \(O(Kd^2+d^3)\), while the trained actor is evaluated once per outer decision epoch with fixed cost \(C_{\pi}\). Thus, for \(M\) detector samples and \(T\) decision epochs, the total online cost is \(O(M(Kd^2+d^3)+TC_{\pi})\). In the stepper-motor implementation, \(d=1\) and \(K=2000\) are fixed, so the detector update has constant per-sample cost and total runtime scales linearly with the number of acquired samples. Memory is \(O(M)\) when logs are stored and can be reduced to \(O(1)\) with a sliding-window implementation.

4.5.2 Performance Evaluation↩︎

Figures 11 (b)–11 (e) show that the trained RL policy successfully transfers to the physical testbed. Under the replay attack starting at decision epoch \(7\) (iteration \(\sim600\)), the physical trajectory in Fig. 11 (c) diverges from the spoofed measurements, while the state estimate \(\widehat{\mathbf{y}}_t\) tracks the replayed input. The discrepancy between the spoofed measurement and state estimates triggers the detector belief \(d_t\) in Fig. 11 (d), which reaches \(1\) within five samples. Figure 11 (e) shows DynaMark adapting \(U_t\): under nominal motion, it peaks near \(\rm 0.009 \;mV^2\) during acceleration and tapers on descent; after belief updating to 1, DynaMark adjusts \(U_t\) to maintain a strong residual gap. These results confirm that DynaMark is deployable on the physical implementation while successfully distinguishing replayed measurements from the actual motor response in real time.

4.5.3 Comparative Analysis↩︎

To instantiate the LTI-based baselines for the stepper-motor DE, we form a single LTI surrogate despite the mode-switching behavior shown in Fig. 11 (a). Specifically, we estimate an ARX(1,1) model from watermark-free data over the full motion cycle, obtaining \(A=1.0\), \(B=0.00372\), and \(Q=3.1911{\times}10^{-5}\). Using this surrogate, we compute offline static watermark variances by solving the covariance-design optimization in [8] under loss budgets \(\delta\in\{100\%,50\%,20\%,10\%,5\%\}\) relative to the no-watermark LQG cost \(J_0\) with \(X=1.0\), and we compute the online baseline by implementing the adaptive method in [8] and reporting its converged variance. Table ¿tbl:tab:opt95u? lists the resulting \(U_\star\). Each baseline is then evaluated on the DE under the same nominal and replay settings as in Section 4.5.1.

1. For the offline LTI-designed static watermark, tightening \(\delta\) reduces watermark energy and control-performance degradation, but the wide inter-quartile ranges indicate sensitivity to the time-varying motor regime. The online LTI baseline operates with smaller nominal cost across all budgets, but DynaMark maintains consistently low watermark energy and nominal performance degradation without relying on a static covariance tuned to an LTI surrogate.

2. The post-onset \(\mathrm{AF}\) decreases as \(\delta\) tightens for the offline baseline, showing the expected trade-off between nominal cost and detection performance. The online LTI baseline yields near-zero \(\mathrm{AF}\) across budgets, suggesting that its converged watermark intensities are too weak for reliable replay detection on the time-varying motor trajectory. DynaMark achieves substantially higher \(\mathrm{AF}\) with no dependence on \(\delta\), indicating stronger attack sensitivity while remaining performance-aware in nominal operation.

5 Conclusion↩︎

This paper developed DynaMark, an RL-based DWM framework to detect replay attacks in industrial MTCs. By casting adaptive watermark design as an MDP, DynaMark learns an adaptive covariance policy from measured data and detector feedback without supplying an analytic plant transition model during policy optimization. The detector and belief-update modules are model-augmented: under Gaussian residuals, theoretical results characterize how input watermarking preserves the nominal \(\chi^{2}\) detector structure and yields a tractable generalized \(\chi^{2}\) statistic under replay attack, enabling a Bayesian belief update used as a compact detection-confidence state. In addition, under locally contractive nominal closed-loop operation and bounded watermark covariance, we established a mean-square boundedness guarantee for the tracking error, providing a principled safety envelope for learned watermark policies.

Experiments on the Siemens Sinumerik 828D DE, the nonlinear MSD benchmark, and a physical stepper-motor testbed show that DynaMark learns an adaptive covariance policy \(U_t\) from measurements and detector beliefs, improving the detection–performance trade-off without relying on a static watermark intensity. In the DE and MSD studies, DynaMark outperforms both static-variance baselines and a non-RL belief-adaptive baseline, showing that reward-optimized covariance adaptation is more reliable under replay than direct belief-based scaling. The MSD case further demonstrates applicability under time-varying and non-Gaussian conditions by replacing the Gaussian closed-form belief computation with offline Monte-Carlo calibration, without invoking the closed-loop boundedness theorem. On the physical stepper-motor, DynaMark transfers from DE to hardware and maintains post-onset alarm activity under mode-switching dynamics, whereas LTI-surrogate baselines can become either too costly or too weak for reliable detection. Future work will integrate DynaMark with online watermark recovery to synthesize authenticated, energy-efficient recovery inputs after detection.

Acknowledgments↩︎

This work was supported by the U.S. National Science Foundation under Grant 2622141.

5.1 Proof of Lemma 1↩︎

By the replay dynamics in Eq. 8 , the detector forms \(\widehat{\mathbf{y}}'_{t+1}=\mathcal{G}_t(\mathbf{y}'_t,\mathbf{u}^\phi_t)\). Therefore, \(\mathbf{r}_{t+1}^A = \mathbf{e}'_t + \left(\mathcal{G}_t(\mathbf{y}'_t,\mathbf{u}^{\phi'}_t)-\mathcal{G}_t(\mathbf{y}'_t,\mathbf{u}^\phi_t)\right)\). Using Eq. 9 in the input argument around \((\widehat{\mathbf{y}}_t,\Bar{\mathbf{u}}_t)\) yields \(\mathcal{G}_t(\mathbf{y}'_t,\mathbf{u}^{\phi'}_t)-\mathcal{G}_t(\mathbf{y}'_t,\mathbf{u}^\phi_t)\approx H_t(\mathbf{u}^{\phi'}_t-\mathbf{u}^\phi_t)\). Substituting Eq. 6 gives \(\mathbf{r}_{t+1}^A\approx \mathbf{e}'_t + H_t(\phi'_t-\phi_t)\). By distribution of \(\mathbf{e}'_t\), \(\phi_t\), and \(\phi'_t\), and their independence, the affine Gaussian form implies \(\mathbf{r}_{t+1}^A\) is (approximately) Gaussian with zero mean and covariance \(Q_t + H_t(U'_t+U_t)H_t^\top\).0◻

5.2 Proof of Theorem 1↩︎

Under a replay attack, Lemma 1 shows \(\mathbf{r}_{t}^A\sim\mathcal{N}(0, \mathcal{S}_{t})\). We prove the theorem for a generic case where \(\mathbf{r}_{t}^A\sim\mathcal{N}(\mathbf{m}_{t}, \mathcal{S}_{t})\). Using standard normal random variables, we can rewrite \(\mathbf{r}_{t}=\mathbf{m}_{t} + \mathcal{S}_{t}^{1/2} \mathbf{Z}\) where \(\mathbf{Z}\sim\mathcal{N}(0,I)\). Then, the test statistics under attack follow \[\begin{align} g_{t}^A &= \left(\mathbf{Z} + \mathcal{S}_{t}^{-1/2} \mathbf{m}_{t}\right)^\top \mathcal{S}_{t}^{1/2} Q_{t-1}^{-1} \mathcal{S}_{t}^{1/2} \left(\mathbf{Z} + \mathcal{S}_{t}^{-1/2} \mathbf{m}_{t}\right) . \end{align}\] Applying the spectral theorem, \(\mathcal{S}_{t}^{1/2} Q_{t-1}^{-1} \mathcal{S}_{t}^{1/2} = P_{t}^\top \Lambda_{t}P_{t}\) where \(P_{t}\) is an orthogonal matrix, i.e., \(P_{t}^\top P_{t}= P_{t}P_{t}^\top =I\); and, \(\Lambda_{t}\) is a diagonal matrix with positive diagonal elements. Let \(\mathbf{Y}_{t} = P_{t}\mathbf{Z}\) and \(\mathbf{b}_{t}=P_{t}\mathcal{S}_{t}^{-1/2} \mathbf{m}_{t}\). Since \(\mathbf{Z}\) is a multivariate standard normal distribution, \(\mathbf{Y}_{t}\) also follows the same distribution. Then, \[\begin{align} g_{t}^A =& (\mathbf{Y}_{t} + \mathbf{b}_{t})^\top \Lambda_{t} (\mathbf{Y}_{t} + \mathbf{b}_{t}) \\ =& \sum_{i=1}^{n} \Lambda_{t}(i) \left(\mathbf{Y}_{t}(i) + \mathbf{b}_{t}(i)\right)^{2}, \end{align}\] that implies, test statistics under replay follow the generalized \(\chi^{2}\) distribution (weighted sum of noncentral chi-square variables) \(g_{t}^A\sim\tilde{\chi}(\omega_{t},\kappa, \lambda_{t}, s, m)\) where \(s=0\), \(m=0\), \(\kappa =\Vec{\pmb{1}}\), \(\omega_{t} =\left(\Lambda_{t}(1),\dots,\Lambda_{t}(n)\right)^\top\), and \(\lambda_{t}= \mathbf{b}_{t}\). Moreover, by Lemma 1, \(\mathbf{m}_{t}=\mathbf{0}\) (hence \(\mathbf{b}_{t}=\mathbf{0}\)), so \(g_{t}^A\) is a central generalized \(\chi^2\) (a weighted sum of independent \(\chi^2\) variables).0◻

5.3 Proof of Theorem 2↩︎

By the definition of the random variables \(\sigma\) and \(I_t\), the detector’s belief at time \(t\) is defined as \(d_t = \Pr(\sigma=1|I_{1:t})\). Observing \(\mathbf{y}_{t}\) and \(I_t\), and applying the Bayesian rule, \[\begin{align} p(\sigma | I_{1:t}) &= p(\sigma|I_{1:t-1}) \frac{p(I_{t}|I_{1:t-1},\sigma)}{p(I_{t}|I_{1:t-1})}. \end{align}\] Conditioning on \(\sigma\), the conditional evidence is computed as \[\begin{align} p(I_{t}|I_{1:t-1}) =& \sum_{k\in\{0,1\}} p(I_{t},\sigma=k|I_{1:t-1}) \\ =& \sum_{k\in\{0,1\}} p(I_{t}|I_{1:t-1},\sigma=k) \Pr(\sigma=k|I_{1:t-1}). \end{align}\] Under our assumptions, \(\mathbf{e}_t \perp \mathbf{e}_s\) for \(t\ne s\). When \(\sigma=0\), we have \(\mathbf{r}_{t+1}=\mathbf{e}_t\) that implies \(\mathbf{r}_t \perp \mathbf{r}_s\) for \(t\ne s\) given \(\sigma=0\); this implies \(g_t\perp g_s\), and, accordingly, \(I_t\perp I_s\) for \(t\ne s\) given \(\sigma=0\). Therefore, conditional independence under the no-attack hypothesis, the conditional PDF simplifies to \(p(I_{t}|I_{1:t-1},\sigma=0) = p(I_{t}|\sigma=0).\) Therefore, the probability of raising an alarm depend only on the Type-I error, i.e., \(I_{t}|\sigma=0\sim \text{Bernoulli}(\alpha)\).

Define the random variable \(\tau\in\mathbb{R}_{>0}\) as the attack onset. Condition on whether the attack has started by time \(t\): \[\begin{align} p(I_{t} &\mid I_{1:t-1},\sigma=1) = \notag\\ &p(I_{t}\mid I_{1:t-1},\sigma=1, \tau > t) \;p(\tau > t \mid I_{1:t-1},\sigma=1) \notag\\ &+ p(I_{t}\mid I_{1:t-1},\sigma=1, \tau \leq t) \;p(\tau \leq t \mid I_{1:t-1},\sigma=1) \label{eq:prob1} \end{align}\tag{12}\] For \(\tau > t\), since attack has not started yet, the residuals follow the nominal ones. Therefore, by \(\mathbf{r}_t \perp \mathbf{r}_s\), and, accordingly, \(I_t\perp I_s\) for \(t\ne s\) given \(\sigma=1\) and \(\tau > t\), it implies \(p(I_{t}\mid I_{1:t-1},\sigma=1, \tau > t) = p(I_{t}|\sigma=0)\). When \(\tau \leq t\), the attack has started and residuals follow the Gaussian distribution in Lemma 1. Therefore, by Theorem 1, definition of \(\mathbf{r}_t^A\), \(\mathbf{e}_t \perp \mathbf{e}_s\) and \(\phi_t \perp \phi_s\) for \(t \ne s\), and \(\mathbf{e}_t \perp \phi_t\) for all \(t\), \(p(I_{t}\mid I_{1:t-1},\sigma=1, \tau \leq t) = p(I_t \mid \sigma=1)\). Therefore, depending on \(\tau\), probability of raising an alarm depends on both Type-I and Type-II errors, i.e., \[I_{t}|\sigma=1,\tau=k \sim \begin{cases} \text{Ber}(\alpha) & \text{if } t < k\\ \text{Ber}(1 - \beta_{t}) & \text{if } t \geq k \end{cases}.\] Substituting these results in Eq. 12 yields \[\begin{align} p(I_{t} &\mid I_{1:t-1},\sigma=1) = \notag\\ &\alpha^{I_{t}} (1 - \alpha)^{1-I_{t}} \;p(\tau > t \mid I_{1:t-1},\sigma=1) \notag\\ &+ (1 - \beta_{t})^{I_{t}} \beta_{t}^{1-I_{t}} \;p(\tau \leq t \mid I_{1:t-1},\sigma=1), \end{align}\] that yields the conditional evidence \(p(I_{t}|I_{1:t-1})\) and \(d_t\). Before the attack onset, the alarm sequence \(I_{1:t}\) contains only false alarms; after the attack onset, it no longer affects the onset time. Thus, the attack onset time is conditionally independent of the alarm sequence given \(\sigma = 1\), i.e., \(p(\tau \leq t \mid I_{1:t-1}, \sigma = 1) = p(\tau \leq t \mid \sigma = 1).\)

Assuming a prior distribution on the random variable \(\sigma\), which by its definition is a Bernoulli distribution with some parameter \(q\), and the random variable \(\tau|\sigma=1\), the detector’s belief updates at each time using the relation in Eq. ?? .0◻

5.4 Proof of Lemma 2↩︎

Note that \(\alpha=\Pr(I_t=1|\sigma=0)\) and \(\beta_t=\Pr(I_t=0|\sigma=1)\). It is assumed that the Type-I error is time-invariant and can be controlled by carefully selecting the test statistic threshold, \(\tilde{g}\). However, an operator could choose to vary \(\tilde{g}\) over time, which would make the Type-I error time-dependent. The Type-II error, on the other hand, can vary depending on the occurrence of a replay attack. Let the random variable \(\tau\in\mathbb{R}_{>0}\) be the attack onset. Given \(\sigma=1\), if \(\tau \leq t\), then by Theorem 1, the probability of not raising an alarm depends on \(\tilde{\chi}(\omega_{t},\kappa, \lambda_{t}, s, m)\). If \(\tau > t\), the probability of not raising an alarm is based on the \(\chi^2\) distribution. By law of total probability, this relationship is given by \[\begin{align} \beta_t =& \Pr(I_t=0|\sigma=1) \\ =&\Pr(I_t=0, \tau\leq t|\sigma=1) + \Pr(I_t=0, \tau > t|\sigma=1) \\ =& \Pr(g_{t} < \tilde{g} | \sigma = 1, \tau\leq t) \Pr(\tau\leq t | \sigma = 1) \\ &+ \Pr(g_{t} < \tilde{g}) \Pr(\tau > t | \sigma = 1) \\ =& \mathcal{H}_{g_t^A}(\tilde{g}) \Pr(\tau\leq t | \sigma = 1) + \mathcal{H}_{g_t}(\tilde{g}) \Pr(\tau > t | \sigma = 1), \end{align}\] where \(\mathcal{H}_{g_t^A}\) is the CDF of the (central) generalized \(\chi^2\) distribution, and \(\mathcal{H}_{g_t}\) is the CDF of the \(\chi^2\) distribution.0◻

5.5 Estimating \(\beta_t\) under non-Gaussian noise↩︎

Although a Gaussian distribution was assumed, it remains feasible to obtain an estimate of \(\beta_t\). One possible approach is to conduct numerical simulations, provided that a sufficiently accurate digital twin of the system is available. Procedure 13 estimates the lookup table \(\{\hat{\beta}_t(U)\}\) offline for a system under non-Gaussian noise. The detector threshold \(\tilde{g}\) is first calibrated from nominal rollouts to achieve the target false-alarm level \(\alpha\). Then, for each candidate watermark variance \(U\) and replay onset \(\tau\), replay-attack simulations are performed to estimate the probability of no alarm at each time step. Averaging over the chosen onset-time distribution yields \(\hat{\beta}_t(U)\), which replaces \(\beta_t\) in the online belief update.

During online learning, if the selected action \(U_t\) coincides with a grid point in \(\mathcal{U}\), then \(\hat{\beta}_t(U_t)\) is read directly from the lookup table. Otherwise, if \(U^{(k)} \leq U_t \leq U^{(k+1)}\), we use piecewise linear interpolation: \[\begin{align} \hat{\beta}_t(U_t) =& \frac{U^{(k+1)}-U_t}{U^{(k+1)}-U^{(k)}} \, \hat{\beta}_t\!\left(U^{(k)}\right) \\ &+ \frac{U_t-U^{(k)}}{U^{(k+1)}-U^{(k)}} \, \hat{\beta}_t\!\left(U^{(k+1)}\right). \end{align}\] The interpolated value \(\hat{\beta}_t(U_t)\) is then substituted into \(\kappa_{1,t}\) in Theorem 2.

If the watermark action is scalar, or the covariance matrix is parameterized by a scalar design variable, then \(\hat{\beta}_t(U_t)\) can be computed by piecewise linear interpolation over the calibrated action grid. For matrix-valued actions, the same idea applies after replacing \(U_t\) with its low-dimensional parameterization; otherwise, multidimensional interpolation would be required.

5.6 Proof of Theorem 3↩︎

Under the assumptions in § 3.2, let \(\mathbf{y}^\star_t\) be the nominal trajectory and \(\mathbf{u}^\star_t = f_c(\mathbf{h}^\star_t;\eta)\) the ideal control input that keeps the system on this path, i.e., \(\mathbf{y}^\star_{t+1} = \mathcal{G}_t(\mathbf{y}^\star_t, \mathbf{u}^\star_t) + \mathbf{e}^\star_t\). Define the tracking error \(\xi_t=\mathbf{y}_t - \mathbf{y}^\star_t\), and, the control deviation \(\Delta\mathbf{u}_t = \mathbf{u}_t - \mathbf{u}^\star_t\). In a neighborhood of \((\mathbf{y}^\star_t,\mathbf{u}^\star_t)\), a first-order Taylor expansion gives \[\label{eq:taylor95g} \mathcal{G}_t(\mathbf{y},\mathbf{u})\approx \mathcal{G}_t(\mathbf{y}^\star_t,\mathbf{u}^\star_t) + G_t^\star(\mathbf{y}-\mathbf{y}^\star_t) + H_t^\star(\mathbf{u}-\mathbf{u}^\star_t),\tag{13}\] where higher-order terms are negligible in the operating region of interest, \(G_t^\star:=\nabla_y \mathcal{G}_t(\mathbf{y}^\star_t,\mathbf{u}^\star_t)\) and \(H_t^\star:=\nabla_u \mathcal{G}_t(\mathbf{y}^\star_t,\mathbf{u}^\star_t)\). By differentiability of \(\mathcal{G}_t\) on a compact neighborhood of \((\mathbf{y}^\star_t,\mathbf{u}^\star_t)\) for all \(t\), there exists \(\Bar{H}>0\) such that \(\sup_{t\ge0} \|H_t^\star\|^2\leq\Bar{H}\), for all \((\mathbf{y},\mathbf{u})\) in this neighborhood. By differentiability of the controller in a neighborhood of \(\mathbf{h}_t^\star\), we obtain the linear approximation \(f_c(\mathbf{h};\eta) \approx f_c(\mathbf{h}_t^\star;\eta) + K_t(\mathbf{h}- \mathbf{h}_t^\star)\), where \(K_t = [K_{t,0},\dots,K_{t,w}] \in \mathbb{R}^{c\times n(w+1)}\) and \(K_{t,i} = \nabla_{\mathbf{y}_{t-i}} f_c(\mathbf{h}_t^\star;\eta)\) for \(i=0,\dots,w\). This yields \(\mathbf{h}- \mathbf{h}_t^\star = \left[\xi_t, \xi_{t-1}, \dots, \xi_{t-w}\right]^\top\) and \(\Delta\mathbf{u}_t \approx \sum_{i=0}^{w} K_{t,i} \xi_{t-i}\). Since the nominal histories lie in a compact set, then \(\sup_{t>0}\|K_{t,i} \|^2 \leq \Bar{K}_i\) for all \(i=0,\dots,w\).

Then, using Eq. 13 , locally around \((\mathbf{y}^\star_t,\mathbf{u}^\star_t)\), we have \(\xi_{t+1} \approx A_t \xi_t + d_t,\) where \(A_t = G_t^\star + H_t^\star K_{t,0}\), \(d_t = H_t^\star \sum_{i=1}^{w} K_{t,i} \xi_{t-i} + H_t^\star \phi_t + \varrho_t\), and \(\varrho_t=\mathbf{e}_t - \mathbf{e}^\star_t + \text{unmodeled noises}\). For any \(\epsilon>0\), Young’s inequality gives: \[\label{eq:ineq1} \|\xi_{t+1}\|^2 \leq (1+\epsilon) \|A_t \xi_t\|^2 + \left(1+\epsilon^{-1}\right) \|d_t\|^2.\tag{14}\]

By Cauchy–Schwarz (CS) and boundedness of Jacobians, \(\|d_t\|^2 \leq 3 \left(w\Bar{H} \sum_{i=1}^{w} \Bar{K}_i \|\xi_{t-i}\|^2 + \Bar{H} \|\phi_t\|^2 + \|\varrho_t\|^2\right)\) that depends on \(\{\xi_\tau\}_{\tau=t-w}^{t}\). Define \[\label{eq:lyp} V_t = \sum\nolimits_{i=0}^{w} \rho^i \;\|\xi_{t-i}\|^2 , \qquad \rho\in(0,1).\tag{15}\] By Eq. 15 , \(\|\xi_{t-i}\|^2 \leq \rho^{-i} V_t\). Hence, \(\sum_{i=1}^{w} \Bar{K}_i\|\xi_{t-i}\|^2 \leq V_t \sum_{i=1}^{w} \Bar{K}_i \rho^{-i}.\) Moreover, with \(\Bar{K}=\max_{i=1,\dots,w}\Bar{K}_i\) and \(\rho\in(0,1)\), \(\sum_{i=1}^{w} \Bar{K}_i \rho^{-i} \leq w\Bar{K} \rho^{-w}\). Combining together, \[\label{eq:ineq3} \|d_t\|^2 \leq 3 \left(w^2 \Bar{H} \Bar{K} \rho^{-w} V_t + \Bar{H} \|\phi_t\|^2 + \|\varrho_t\|^2\right).\tag{16}\]

By CS, \(\|A_t \xi_t\|^2 \leq \|A_t \|^2 \;\|\xi_t\|^2\), and by Eq. 15 , \(\|\xi_t\|^2 \leq V_t\). Then, \(\|A_t \xi_t\|^2 \leq \|A_t \|^2 V_t\). Together with 14 and 16 \[\begin{align} \label{eq:ineq4} \|\xi_{t+1}\|^2 \leq& \left((1+\epsilon) \|A_t \|^2 + 3\left(1+\epsilon^{-1}\right)w^2 \rho^{-w} \Bar{H} \Bar{K} \right)V_t \notag \\ &+ 3\left(1+\epsilon^{-1}\right) \left(\Bar{H} \|\phi_t\|^2 + \|\varrho_t\|^2\right). \end{align}\tag{17}\]

By Eq. 15 , \[\begin{align} V_{t+1} =\|\xi_{t+1}\|^2 + \rho \sum_{i=0}^{w-1} \rho^{i} \;\|\xi_{t-i}\|^2\leq \|\xi_{t+1}\|^2 + \rho V_t.\label{eq:ineqv} \end{align}\tag{18}\] Adding \(\rho V_t\) to both sides of 17 , using 18 , and taking the expectation yield \[\begin{align} \mathop{\mathrm{\mathbb{E}}}[V_{t+1}] \leq& \; \left((1+\epsilon) \|A_t \|^2 + \rho + 3\left(1+\epsilon^{-1}\right)w^2 \rho^{-w} \Bar{H} \Bar{K} \right) \mathop{\mathrm{\mathbb{E}}}[V_t] \\ &+ 3\left(1+\epsilon^{-1}\right) \left(\Bar{H} \mathop{\mathrm{\mathbb{E}}}\left[\|\phi_t\|^2\right] + \mathop{\mathrm{\mathbb{E}}}\left[\|\varrho_t\|^2\right]\right). \end{align}\] As \(\phi_t\sim\mathcal{N}(0,U_t)\) and \(U_t\in\mathcal{A}\), by CS and  10 , we have that \(\mathop{\mathrm{\mathbb{E}}}\left[\|\phi_t\|^2\right] = \mathop{\mathrm{tr}}(U_t) = \left\langle U_t,I\right\rangle \leq \|U_t\|_F \|I\|_F \leq U_{\max}\sqrt{c}\). Also, assume \(\mathop{\mathrm{\mathbb{E}}}\left[\|\varrho_t\|^2\right] \leq \sigma^2_\varrho\). Altogether, \[\begin{align} \mathop{\mathrm{\mathbb{E}}}[V_{t+1}] \leq& \left((1+\epsilon) \|A_t \|^2 + \rho + 3\left(1+\epsilon^{-1}\right)w^2 \rho^{-w} \Bar{H} \Bar{K} \right) \mathop{\mathrm{\mathbb{E}}}[V_t] \\ &+ 3\left(1+\epsilon^{-1}\right) \left(\Bar{H} U_{\max}\sqrt{c} + \sigma^2_\varrho\right). \end{align}\] Let \(\Bar{A}=\sup_{t\ge0} \|A_t \|\). Then, \[\label{eq:ineq6} \mathop{\mathrm{\mathbb{E}}}[V_{t+1}] \le \nu(\epsilon,\rho)\mathop{\mathrm{\mathbb{E}}}[V_t] + \vartheta(\epsilon, U_{\max}),\tag{19}\] where \[\begin{align} \nu(\epsilon,\rho) &= (1+\epsilon) \Bar{A}^2 + \rho + 3\left(1+\epsilon^{-1}\right)w^2 \rho^{-w} \Bar{H} \Bar{K}, \\ \vartheta(\epsilon, U_{\max}) &= 3\left(1+\epsilon^{-1}\right) \left(\Bar{H} U_{\max}\sqrt{c} + \sigma^2_\varrho\right). \end{align}\] To ensure 19 is a contraction, it must \(\nu(\epsilon,\rho)<1\).

Assume the nominal closed loop is locally contractive, i.e., \(\Bar{A}<1\). Then, there exists \(\epsilon>0\) such that \((1+\epsilon)\Bar{A}^2<1\). Define the residual margin \(\delta_\epsilon = 1-(1+\epsilon)\Bar{A}^2>0\) and let \[\label{eq:c95epsilon} c_\epsilon=3(1+\epsilon^{-1})w^2 \Bar{H}\Bar{K}.\tag{20}\] A sufficient condition for \(\nu(\epsilon,\rho)<1\) is the existence of \(\rho\in(0,1)\) with \(g_\epsilon(\rho)=\rho + c_\epsilon \rho^{-w} < \delta_\epsilon\). \(g_\epsilon\) is strictly convex on \((0,1)\) and has a unique minimizer \(\rho_\epsilon^\star\) (if \(\rho_\epsilon^\star\in(0,1)\)). Thus, it suffices that some \(\rho_\epsilon^\star\in\arg\min_{\rho\in(0,1)} g_\epsilon(\rho)\) satisfies \(g_\epsilon(\rho_\epsilon^\star)<\delta_\epsilon\). Differentiating \(g_\epsilon\) and setting the derivative to zero gives \[\label{eq:rho95g95star} \rho_\epsilon^\star = (w c_\epsilon)^{\frac{1}{w+1}}, \qquad g_\epsilon(\rho_\epsilon^\star) = (1+w^{-1})\rho_\epsilon^\star.\tag{21}\] Note that, by definition of \(\delta_\epsilon\), \(g_\epsilon(\rho_\epsilon^\star)<\delta_\epsilon\) implies \(\rho_\epsilon^\star\in(0,1)\). Therefore, we have the following conditions for \(\nu(\epsilon,\rho)<1\):

• Existence of \(\epsilon>0\): A necessary and sufficient condition is \((1+\epsilon)\Bar{A}^2 < 1\) for some \(\epsilon>0\), which holds iff \(\Bar{A}<1\) (e.g., any \(\epsilon\in(0,\Bar{A}^{-2}-1)\) is admissible). Without this, \(\nu(\epsilon,\rho)<1\) is not guaranteed.

• Existence of \(\rho\in(0,1)\) for a given \(\epsilon>0\): It suffices that \(g_\epsilon(\rho_\epsilon^\star)<\delta_\epsilon\).

For the latter condition, let \(\epsilon_\alpha=\alpha(\Bar{A}^{-2} - 1)\) be admissible, with \(\alpha \in (0,1)\) and \(\Bar{A} < 1\). Substituting \(\epsilon_\alpha\) into Eq. 20 yields the corresponding coefficient \(c_{\epsilon_\alpha}\). In turn, substituting \(c_{\epsilon_\alpha}\) into the expression for \(\rho_{\epsilon}^\star\) in Eq. 21 leads to \[\rho_{\epsilon_\alpha}^\star = \left(3 w^3 \Bar{H} \Bar{K} \,\frac{\alpha + (1-\alpha)\Bar{A}^2}{\alpha(1 - \Bar{A}^2)}\right)^{\frac{1}{w+1}}.\] By \(\alpha + (1-\alpha)\Bar{A}^2 < 1\), it follows that \[\rho_{\epsilon_\alpha}^\star \;<\; \left(\frac{3w^3 \Bar{H} \Bar{K}}{\alpha\bigl(1 - \Bar{A}^2\bigr)}\right)^{\frac{1}{w+1}}.\] Therefore, by the expression for \(g_\epsilon(\rho_\epsilon^\star)\) in Eq. 21 and with some \(\alpha\in(0,1)\), if \[\label{eq:ineq7} \left(1+w^{-1}\right)\left(\frac{3w^3 \Bar{H} \Bar{K}}{\alpha(1 - \Bar{A}^2)}\right)^{\frac{1}{w+1}} < \delta_{\epsilon_\alpha},\tag{22}\] then defined \(\epsilon\) satisfies \(g_{\epsilon}(\rho_{\epsilon_{\alpha}}^\star) < \delta_{\epsilon_\alpha}\) that implies \(\rho\) exists, and \(\rho=\rho^\star_\epsilon\). Substituting \(\epsilon_\alpha\) in the expression for \(\delta_{\epsilon_\alpha}\) yields \(\delta_{\epsilon_\alpha}=(1-\alpha)(1-\Bar{A}^2)\). In turn, substituting \(\delta_{\epsilon_\alpha}\) in 22 , the inequality 22 holds if \[\tag{23} \begin{align} &\Bar{A} < \sqrt{1 - \left(\frac{3w^3 \left(1+w^{-1}\right)^{w+1} \Bar{H} \Bar{K}}{\alpha(1-\alpha)^{w+1}}\right)^{\frac{1}{w+2}}}, \tag{24} \\ &3w^3 \left(1+w^{-1}\right)^{w+1} \Bar{H} \Bar{K} < \alpha(1-\alpha)^{w+1}. \end{align}\] Define \(f_\epsilon(\alpha)=\alpha(1-\alpha)^{w+1}\). The largest \(\alpha\) for which 23 holds is the maximizer of \(f_\epsilon\). \(f_\epsilon\) is concave on \((0,1)\), and differentiating and setting the derivative to zero gives \(\alpha^\star=\frac{1}{w+2}\), with \(f_\epsilon(\alpha^\star)=\frac{1}{w+2}\left(\frac{w+1}{w+2}\right)^{w+1}\). Substituting \(f_\epsilon(\alpha^\star)\) in 23 , the necessary and sufficient conditions for the existence of \(\epsilon>0\) and \(\rho\in(0,1)\) are (since 24 is stricter than \(\bar{A}<1\)): \[\label{eq:ineq11} \begin{align} \Bar{A} &< \sqrt{1 - \left(\frac{3 (w+2)^{w+2} \Bar{H} \Bar{K}}{w^{w-2}}\right)^{\frac{1}{w+2}}}, \\ \Bar{H} \Bar{K} &< \frac{w^{w-2}}{3(w+2)^{w+2}}. \end{align}\tag{25}\]

Thus, under the conditions in 25 with \(\nu(\epsilon,\rho)<1\), 19 yields a linear recursion: \(\mathop{\mathrm{\mathbb{E}}}[V_{t}] \le \nu(\epsilon,\rho)^t V_0 + \left[\vartheta(\epsilon,U_{\max})/(1-\nu(\epsilon,\rho))\right]\), and, by \(\|\xi_{t}\|^2 \leq V_t\), it yields \(\mathop{\mathrm{\mathbb{E}}}\|\xi_{t}\|^2 \le \nu(\epsilon,\rho)^t V_0 + \left[\vartheta(\epsilon,U_{\max})/(1-\nu(\epsilon,\rho))\right]\). Therefore, \[\limsup_{t\to\infty}\mathop{\mathrm{\mathbb{E}}}\|\xi_{t}\|^2 \leq \frac{\vartheta(\epsilon,U_{\max})}{1-\nu(\epsilon,\rho)}. \tag*{\qed}\]

5.6.1 Theorem 3 Verification↩︎

Given the controller history window \(w\), verify Theorem 3 as follows:

1. Compute \(\Bar{H}\), \(\Bar{K}\), and \(\Bar{A}\) from the local closed-loop model.

2. Verify \(\mathop{\mathrm{\mathbb{E}}}\|\varrho_t\|^2\le \sigma_\varrho^2<\infty\).

3. Check conditions \((C_1)\) and \((C_2)\). If either condition fails, the existence of admissible \((\epsilon,\rho)\) is not guaranteed by Theorem 3.

4. Choose \(\epsilon\in(0,\Bar{A}^{-2}-1)\) and compute \(\delta_\epsilon\) and \(c_\epsilon\).

5. If \(c_\epsilon>0\), compute \(\rho_\epsilon^\star\) and set \(\rho=\rho_\epsilon^\star\). Then verify \(g_\epsilon(\rho_\epsilon^\star)<\delta_\epsilon\) and \(\nu(\epsilon,\rho)<1\).

6. If \(c_\epsilon=0\), choose any \(\rho\in(0,\delta_\epsilon)\) that gives \(\nu(\epsilon,\rho)<1\).

7. Specify \(U_{\max}\), compute \(\vartheta(\epsilon,U_{\max})\), and report the resulting mean-square bound in ?? .

5.7 Implementation Specifications↩︎

To deploy the DynaMark framework, we implemented a DDPG-CDQ agent using a PyTorch-based actor-critic setup. Table ¿tbl:tab:ddpg95specs? lists the DDPG-CDQ hyperparameters for our experiments. Implementation details, training hyperparameters, and code to reproduce the experiments are available at: [https://github.com/navidaftabi/DynaMark]. Figure 14 reports representative training and evaluation returns for the DDPG-CDQ used to instantiate DynaMark in the three case studies. These curves are provided as empirical training diagnostics, not as a theoretical convergence guarantee. other algorithms such as PPO, SAC, or full TD3 could be used to solve the same framework. Training stability was supported by clipped double-Q targets, target networks with Polyak averaging, gradient clipping, layer normalization, replay-buffer sampling, and decaying OU exploration. Final policy behavior under reward, watermark-budget, and prior variations is evaluated in the ablation and sensitivity study, where the reported metrics are aggregated over seeds and replications. Reward weights \((\omega_1, \omega_2, \omega_3)\) encode the practitioner’s preference among actuation overhead, nominal performance, and detection confidence. A practical selection procedure is to first set \(\omega_2\) to satisfy a process-specific tracking/quality tolerance under nominal operation (e.g., ensuring \(\|\mathbf{y}_t^{\star} - \mathbf{y}_t\|_2\) remains within an acceptable band), then increase \(\omega_3\) until replay attacks consistently drive the belief away from \(0.5\) with acceptable delay. Finally, tune \(\omega_1\) to limit watermark energy/wear while preserving the achieved detection behavior. In practice, we found that keeping \(\omega_1+\omega_2+\omega_3=1\) and performing a coarse grid search over a small set of candidates yields stable policies, with the final choice guided by the desired operating point on the security–performance frontier.

References↩︎

---

1. Navid Aftabi and Dan Li are with Industrial & Systems Engineering Department, University of Washington, Seattle, WA 98195 USA (e-mail: aftabi@uw.edu, dli27@uw.edu).↩︎

2. Abhishek Hanchate and Satish Bukkapatnam are with Wm Michael Barnes ’64 Department of Industrial and Systems Engineering, Texas A&M University, College Station, TX 77843 USA. (e-mail: abhishek.hanchate@tamu.edu, satish@tamu.edu).↩︎