Closed-loop Diffusion Control of Complex Physical Systems

## Abstract

Abstract The control problems of complex physical systems have broad applications in science and engineering. Previous studies have shown that generative control methods based on diffusion models offer significant advantages for solving these problems. However, existing generative control approaches face challenges in both performance and efficiency when extended to the closed-loop setting, which is essential for effective control. In this paper, we propose an efficient C losed- L oop Diff usion method for Phy sical systems Con trol (CL-DiffPhyCon). By employing an asynchronous denoising framework for different physical time steps, CL-DiffPhyCon generates control signals conditioned on real-time feedback from the environment with significantly reduced computational cost during sampling. Additionally, the control process could be further accelerated by incorporating fast sampling techniques, such as DDIM. We evaluate CL-DiffPhyCon on two tasks: 1D Burgers’ equation control and 2D incompressible fluid control. The results demonstrate that CL-DiffPhyCon achieves superior control performance with significant improvements in sampling efficiency.

## 1 Introduction

The control problem of complex physical systems is a critical area of study that involves optimizing a sequence of control actions to achieve specific objectives.
It has important applications across a wide range of science and engineering fields, including fluid control (Verma et al., ), plasma control (Degrave et al., ), and particle dynamics control (Reyes Garza et al., ). The challenge in controlling such systems arises from their high-dimensional, highly nonlinear, and stochastic characteristics. Therefore, to achieve good control performance, there is an inherent requirement of *closed-loop* control. Specifically, each control decision should be based on the latest state feedback from the environment, allowing for continuous adaptation of the control inputs in response to any changes. This requirement is particularly necessary for control tasks that involve extra challenges, such as partial observations or indirect control.

Over recent decades, several methods have been developed to address this problem, including classical control methods, recent reinforcement learning approaches, and the latest generative methods. Among them, diffusion models, such as DiffPhyCon (Wei et al., ), have demonstrated competitive performance, often outperforming both classical control and reinforcement learning methods in complex physical systems control. The superiority of diffusion models for general decision-making problems has also been widely demonstrated recently (Ajay et al., ; Janner et al., ).

However, these diffusion control approaches encounter significant challenges in handling the closed-loop control problems, due to their reliance on a synchronous denoising strategy.
The diffusion models start from pure noise to a denoised sample for all physical time steps within the model horizon. Applying a full sampling process at each physical time step can realize the closed-loop control, but incurs high sampling cost.
Moreover, it may disrupt the consistency of the control signals, thereby affecting overall performance (Kaelbling & Lozano-Pérez, ).
On the other hand, conducting a full sampling process every several physical time steps improves sampling efficiency but no longer conforms to the closed-loop requirement,
leading to inferior control decisions due to outdated system states.
Although an online replanning strategy has been proposed recently to determine when to replan adaptively (Zhou et al., ), they do not establish a fully closed-loop framework. In addition, it involves extra computation of likelihood estimation or sampling from scratch, with a dependence on thresholding hyperparameters, which may vary across different tasks and require experiments or a significant amount of tuning to determine.

Figure: Figure 1: Advantages of our CL-DiffPhyCon (right) over previous diffusion control methods (left and middle). By employing an asynchronous denoising framework, our method could achieve closed-loop control and accelerate the sampling process significantly. The notation DiffPhyCon-$h$ means conducting a full sampling process including $T$ denoising steps every $h$ physical time steps.
Refer to caption: x1.png

In this paper, we propose a novel Closed-Loop Diffusion method for Physical systems Control, named as CL-DiffPhyCon. The key idea is to decouple the synchronous denoising within the model horizon, allowing different physical time steps to exhibit different noise levels. In this way, online closed-loop generation of control sequences is naturally realized: the asynchronous diffusion model outputs control signals sequentially with increasing levels of noise along physical time steps, which enables real-time interaction between the environment and the sampled control signal in each horizon without waiting for all the following control signals in the same horizon to be denoised completely. Then, the environmental feedback serves as the initial condition for sampling subsequent control signals, ensuring they are generated based on this reliable state. Our method can also be seen as a seamless replanning approach that leverages fresh observations from the environment with minimal sampling costs. Therefore, compared to existing diffusion-based control methods (Wei et al., ; Zhou et al., ), our approach not only realizes closed-loop control but also achieves significant sampling acceleration. These advantages of our method are illustrated in Figure 1.

In summary, we make the following contributions:
(1) We develop CL-DiffPhyCon, a novel closed-loop diffusion control method of complex physical systems. It speeds up diffusion control methods in a closed-loop framework, which is derived from a theoretical analysis of the target distribution we aim to sample.
Additionally, the control process could be further accelerated by incorporating fast sampling techniques, such as DDIM (Song et al., ).
(2) We evaluate CL-DiffPhyCon on the 1D Burgers’ equation control and 2D incompressible fluid control tasks. The results demonstrate that CL-DiffPhyCon achieves notable control performance with significant sampling acceleration.

## 2 Related Work

Classical control methods like PID (Li et al., ) are known for their high efficiency, steady performance, and good interpretability, but they face significant challenges in both performance and efficiency when applied to control high-dimensional complex physical systems. Recently, imitation learning and reinforcement learning (Pomerleau, ; Zhuang et al., ) has shown good performance on a wide range of physical systems, including drag reduction (Rabault et al., ; Elhawary, ; Feng et al., ; Wang et al., ), heat transfer (Beintema et al., ; Hachem et al., ), and fish swimming (Novati et al., ; Verma et al., ; Feng et al., ).
Another category of supervised learning (SL) methods (Holl et al., ; Hwang et al., ) plan control signals by utilizing backpropagation through a neural surrogate model. In contrast, our approach does not depend on surrogate models; instead, it simultaneously learns the dynamics of physical systems and the corresponding control sequences.
Additionally, physics-informed neural networks (PINNs) (Willard et al., ) have recently been utilized for control (Mowlavi & Nabi, ), but they necessitate explicit PDEs. In contrast, our method is data-driven and capable of tackling a wider range of complex physical system control problems provided that their environment is accessible.

Diffusion models (Ho et al., ) excel at learning high-dimensional distributions and have achieved significant success in image and video generation (Dhariwal & Nichol, ; Ho et al., ). It has also demonstrated remarkable capabilities in addressing scientific and engineering challenges, such as weather forecasting (Price et al., ), robot control (Janner et al., ; Ajay et al., ), and inverse problems in PDEs like designing the initial or boundary conditions (Wu et al., ) or planning external temporal control signals (Wei et al., ). However, they are not suitable to be directly applied to efficient closed-loop control since all physical time steps are denoised in the same schedule.
Some previous work has focused on improving the
adaptability of diffusion generation (Zhou et al., ) on decision-making tasks, but it is not closed-loop and needs extra hyperparameters and computations for the decision of replanning. In contrast, our method is a closed-loop approach with an efficient sampling strategy.
Some recent works also assign varying noise levels to different frames within a model horizon (Wu et al., ; Ruhe et al., ). The key differences between our work and these approaches are twofold: (1) Our method focuses on the closed-loop control task of complex physical systems, whereas their methods are geared towards sequential content generation that does not involve interaction with the external world.
(2) The diffusion models we need to learn are derived from the target distribution we aim to sample from (please refer to Section 4.1 for details), while theirs is heuristic.

## 3 Background

### 3.1 Problem Setup

We consider the following online complex physical systems control problem:

$$ $\min_{\mathbf{w}_{1:N}}\mathcal{J}(\mathbf{u}_{0},\mathbf{w}_{1},\mathbf{u}_{1 },...,\mathbf{w}_{N},\mathbf{u}_{N})\quad\text{s.t.}\quad\mathbf{u}_{\tau+1}=G (\mathbf{u}_{\tau},\mathbf{w}_{\tau+1}).$ (1) $$

Here $\mathbf{u}_{\tau}\in\mathbb{R}^{d_{\mathbf{u}}}$ and $\mathbf{w}_{\tau}\in\mathbb{R}^{d_{\mathbf{w}}}$ are the system state and external control signal at physical time step $\tau$, respectively; $G$ is the environment, which is determined by the physical dynamics, like partially differential equations (PDEs).
$\mathcal{J}(\mathbf{u}_{0},\mathbf{w}_{1},\mathbf{u}_{1},...,\mathbf{w}_{N},
\mathbf{u}_{N})$ denotes the control objective defined over a trajectory of length $N$.
For example, $\mathcal{J}$ can be designed to measure the control performance towards a target state $\mathbf{u}^{*}$ with cost constraints: $\mathcal{J}=\|\mathbf{u}_{N}-\mathbf{u}^{*}\|^{2}+\sum_{\tau=1}^{N}\|\mathbf{w
}_{\tau}\|^{2}$. For convenience, we introduce a variable $\mathbf{z}_{\tau}=[\mathbf{w}_{\tau},\mathbf{u}_{\tau}]$ to represent the concatenation of $\mathbf{w}_{\tau}$ and $\mathbf{u}_{\tau}$.
Additionally, we assume the transition probabilities of the system has the following *Markov property*:

$$ $p(\mathbf{z}_{\tau+1}|\mathbf{u}_{0},\mathbf{z}_{1:\tau})=p(\mathbf{z}_{\tau+1 }|\mathbf{u}_{\tau})=p(\mathbf{w}_{\tau+1}|\mathbf{u}_{\tau})p(\mathbf{u}_{ \tau+1}|\mathbf{u}_{\tau},\mathbf{w}_{\tau+1}).$ (2) $$

In this paper, we focus on *closed-loop* control, which means that the control signal $\mathbf{w}_{\tau+1}$ in each time step is planned based on current system state $\mathbf{u}_{\text{env},\tau}$.
The advantage of closed-loop is obvious: real-time feedback from the environment enable the algorithm to adjust its control decisions based on the evolving state of the environment.

### 3.2 Preliminary: Diffusion Control Models

DiffPhyCon (Wei et al., ) is a recent diffusion generative method to solve the problem Eq. 1 for small $N$ in the open-loop manner.
In this section, we briefly review this framework of its light version. Suppose that we have a training set $\mathcal{D}_{\text{train}}$ containing $M$ trajectories $\{\mathbf{u}_{0}^{(i)},\mathbf{z}_{1:N}^{(i)}\}_{i=1}^{M}$ from the environment.
We define the following forward diffusion SDE (Song et al., ) on the $\mathcal{D}_{\text{train}}$,

$$ ${\rm d}\mathbf{z}_{\tau}(t)=f(t)\mathbf{z}_{\tau}(t){\rm d}t+g(t){\rm d}\omega _{\tau}(t),\quad\tau\in[1:N],\quad t\in[0,T],$ (3) $$

where $\mathbf{z}_{\tau}(0)=\mathbf{z}_{\tau}$, $f(t)$ and $g(t)$ are scalar functions, and $\omega_{\tau}(t)$ is Wiener process (^1^11In $\mathbf{z}_{\tau}(t)$, the subscript $\tau$ denotes physical time step, and $(t)$ in the parentheses indicates SDE step.).
Through Eq. 3, we augment the distribution of $\mathcal{D}_{\text{train}}$, denoted as $p_{\text{train}}(\mathbf{u}_{0},\mathbf{z}_{1:N})=p_{\text{train}}(\mathbf{u}_
{0})p_{\text{train}}(\mathbf{z}_{1:N}|\mathbf{u}_{0})$, to the distribution of the diffusion process $p_{t}(\mathbf{u}_{0},\mathbf{z}_{1:N}(t))=p_{\text{train}}(\mathbf{u}_{0})p_{t
}(\mathbf{z}_{1:N}(t)|\mathbf{u}_{0})$.
We have the terminal conditions, $p_{0}(\mathbf{z}_{1:N}(0)|\mathbf{u}_{0})=p_{\text{train}}(\mathbf{z}_{1:N}(0)
|\mathbf{u}_{0})$ and $p_{T}(\mathbf{z}_{1:N}(T)|\mathbf{u}_{0})\approx\mathcal{N}(\mathbf{z}_{1:N}(T
)|\mathbf{0},\sigma_{T}^{2}\mathbf{I})$.
The reverse-time SDE of Eq. 3 has the following formula:

$$ ${\rm d}\mathbf{z}_{1:N}(t)=[f(t)\mathbf{z}_{1:N}(t)-g(t)^{2}\nabla\log p_{t}( \mathbf{z}_{1:N}(t)|\mathbf{u}_{0})]{\rm d}t+g(t){\rm d}\omega_{1:N}(t),\quad t \in[T,0].$ (4) $$

Once we have the score function $\nabla\log p_{t}(\mathbf{z}_{1:N}(t)|\mathbf{u}_{0})$, we can sample $\mathbf{z}_{1:N}(0)\sim p_{\text{train}}(\mathbf{z}_{1:N}(0)|\mathbf{u}_{0})$ from $\mathbf{z}_{1:N}(T)\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ through Eq. 4. According to DiffPhyCon, the control objective $\mathcal{J}(\mathbf{z}_{1:N}(0))$ is achieved by solving the inverse problem (Song et al., ; Chung et al., ):

$$ $\displaystyle{\rm d}\mathbf{z}_{1:N}(t)=$ $\displaystyle[f(t)\mathbf{z}_{1:N}(t)-g(t)^{2}\nabla_{\mathbf{z}_{1:N}(t)}\log p _{t}(\mathbf{z}_{1:N}(t)|\mathbf{u}_{0})$ (5) $\displaystyle+$ $\displaystyle g(t)^{2}\lambda\cdot\nabla_{\mathbf{z}_{1:N}(t)}\mathcal{J}(\hat {\mathbf{z}}_{1:N}(0))]{\rm d}t+g(t){\rm d}\omega_{1:N}(t).$ $$

Here, $\hat{\mathbf{z}}_{1:N}(0)$ is Tweedie’s estimate, $s(t)^{-1}\mathbf{z}_{1:N}(t)+s(t)\sigma(t)^{2}\nabla_{\mathbf{z}_{1:N}(t)}\log
p
_{t}(\mathbf{z}_{1:N}(t)|\mathbf{u}_{0})$, where $s(t)=\exp\big{(}\int_{0}^{t}f(\eta){\rm d}\eta\big{)}$ and $\sigma(t)=\sqrt{\int_{0}^{t}\frac{g(\eta)^{2}}{s(\eta)^{2}}{\rm d}\eta}$.
In practical implementation, these quantities are specialized by VP SDE (Song et al., ), i.e. DDPM (Ho et al., ), and the diffusion model is implemented by a parameterized denoising network $\bm{\epsilon}_{\phi}$ of horizon $H$, which equals $N$.
Please refer to Appendices A and B for more details.

## 4 Method

Figure: Figure 2: CL-DiffPhyCon for closed-loop control. First, it uses the synchronous diffusion model for initialization. Then, it uses the asynchronous diffusion model for online control iteratively. Sampling of each control signal is based on the latest system state feedback from the environment.
Refer to caption: x2.png

In this section, we detail our method CL-DiffPhyCon.
In Section 4.1, we illustrate our idea and derive the two distributions we need to learn.
To sample from them, we present the synchronous and asynchronous diffusion model in Section 4.2 and Section 4.3, respectively.
In Section 4.4, we introduce closed-loop control (inference) by CL-DiffPhyCon, which is also illustrated in Figure 2.

### 4.1 Asynchronously denoising framework

Recall that DiffPhyCon (Wei et al., ) requires latent variables, e.g., system state and control signal, during the reverse diffusion process of horizon $H$ (which is typically much shorter than $N$) being denoised synchronously. As a result, to sample a control signal, a full denoising process of length $T$ over the whole horizon is performed, introducing a large amount of computation over all the latent variables within this horizon.
To address this issue, we propose CL-DiffPhyCon, an *asynchronous denoising process* scheme such that the latent variables of the early physical time step are denoised in advance of the latter ones.
At each physical time step, the sampled control signal is input to the environment and the output state serves as the initial condition for the following denoising process.
Thus, the control signal is planned based on the current system state and closed-loop control is achieved.
Meanwhile, the computational cost between two successive times steps is reduced significantly compared with synchronous sampling.

Formally, we aim to model the joint distribution $p\big{(}\mathbf{z}_{1}(0),\mathbf{z}_{2}(0),\cdots,\mathbf{z}_{N}(0)|\mathbf{u
}_{0}\big{)}$ in a *Markov* complex physical system, where each $\mathbf{z}_{\tau}(0)=[\mathbf{w}_{\tau},\mathbf{u}_{\tau}]$ is a pair of noise-free control signal and system state, and $\mathbf{u}_{0}$ is the initial state.
Denote $\mathbf{z}_{\tau:\tau+H-1}(t)=[\mathbf{z}_{\tau}(t),\mathbf{z}_{\tau+1}(t),
\cdots,\mathbf{z}_{\tau+H-1}(t)]$ as the sequence of hidden variables with the same noise level in the horizontal interval $[\tau:\tau+H-1]$ and $\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)=[\mathbf{z}_{\tau}(t),\mathbf{z}_{\tau+1
}(t+\frac{1}{H}T),\cdots,\mathbf{z}_{\tau+H-1}(t+\frac{H-1}{H}T)]$ as its counterpart with *asynchronous* noise level. Let’s consider the augmented joint distribution

$$ $p\big{(}\mathbf{z}_{1:N}(0),\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T),\cdots, \tilde{\mathbf{z}}_{N+1:N+H}(\frac{1}{H}T)|\mathbf{u}_{0}\big{)}.$ (6) $$

Through sampling from this augmented joint distribution, we can obtain the desired control signals and states sequence $\mathbf{z}_{1:N}(0)\sim p_{\text{train}}\big{(}\mathbf{z}_{1:N}(0)|\mathbf{u}_
{0}\big{)}$. By conditioning on previous variables sequentially, we obtain the following decomposition theorem:

###### Theorem 1 .

Assume that the joint distribution $p\big{(}\mathbf{z}_{1}(0),\mathbf{z}_{2}(0),\cdots,\mathbf{z}_{N}(0)|\mathbf{u
}_{0}\big{)}$ has *Markov property*.
For any $\tau>0$, we assume $\mathbf{z}_{\tau}(T)$ is independently normally distributed with density $\mathcal{N}(\mathbf{z}_{\tau}(T)|\mathbf{0},\sigma_{T}^{2}\mathbf{I})$.
The augmented joint distribution can be decomposed as:

$$ $\displaystyle p\big{(}\mathbf{z}_{1:N}(0),\tilde{\mathbf{z}}_{1:H}(\frac{1}{H} T),\cdots,\tilde{\mathbf{z}}_{N+1:N+H}(\frac{1}{H}T)|\mathbf{u}_{0}\big{)}$ (7) $\displaystyle=$ $\displaystyle p\big{(}\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T)|\mathbf{u}_{0} \big{)}\prod_{\tau=1}^{N}p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf {u}_{\tau-1}(0),\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)} \mathcal{N}(\mathbf{z}_{\tau+H}(T);\mathbf{0},\sigma_{T}^{2}\mathbf{I}).$ $$

The proofs of this theorem and the subsequent propositions are provided in the Appendix C. This theorem implies that to sample from the augmented joint distribution,
we only need to specify two kinds of distributions for ancestral sampling: $p\big{(}\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T)|\mathbf{u}_{0}\big{)}$, and $p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0),\tilde{
\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}$ for $\tau>0$.
Namely, we need one diffusion model to learn the initializing distribution $p\big{(}\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T)|\mathbf{u}_{0}\big{)}$, and another one to learn the transition distribution $p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0),\tilde{
\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}$.

###### Theorem 1 .

### 4.2 Sampling from initializing distribution

To sample from $p\big{(}\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T)|\mathbf{u}_{0}\big{)}$, we train a *synchronous diffusion model* $\bm{\epsilon}_{\phi}$ by the following loss:

$$ $\mathcal{L}_{\text{synch}}=\mathbb{E}_{t,(\mathbf{u}_{0},\mathbf{z}_{1:H}),\bm {\epsilon}}[\|\mathbf{\bm{\epsilon}}-\bm{\epsilon}_{\phi}\big{(}\mathbf{z}_{1: H}(t),\mathbf{u}_{0},t\big{)}\|_{2}^{2}].$ (8) $$

Here, $\mathbf{z}_{1:H}(t)=s(t)\mathbf{z}_{1:H}+s(t)^{2}\sigma(t)^{2}\bm{\epsilon}$, and the expectation is about $t\sim U(0,T)$, $(\mathbf{u}_{0},\mathbf{z}_{1:H})\sim\mathcal{D}_{\text{train}}$ and $\bm{\epsilon}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$.
Since the negative score function minimizes this loss, we get $\nabla_{\mathbf{z}_{1:H}(t)}\log p_{t}\big{(}\mathbf{z}_{1:H}(t)|\mathbf{u}_{0
}\big{)}\approx-\bm{\epsilon}_{\phi^{*}}\big{(}\mathbf{z}_{1:H}(t),\mathbf{u}_
{0},t\big{)}$ with the trained parameters $\phi^{*}$.
After training, we can apply Eq. 4 to sample $\{\mathbf{z}_{1:H}(t)\}$, where $t$ goes from $T$ to $\frac{1}{H}T$, from which we get $\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T)=[\mathbf{z}_{1}(\frac{1}{H}T),\cdots,
\mathbf{z}_{H}(T)]$ following the desired distribution $p\big{(}\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T)|\mathbf{u}_{0}\big{)}$.

### 4.3 Sampling from transition distribution

To sample from $p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0),\tilde{
\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}$, we train an *asynchronous diffusion model* $\bm{\epsilon}_{\theta}$.
Using the notations in Section 3.2, we define an asynchronous forward diffusion SDE:

$$ $\displaystyle{\rm d}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)=$ $\displaystyle[{\rm d}\mathbf{z}_{\tau}(t),{\rm d}\mathbf{z}_{\tau+1}(t+\frac{1 }{H}T),\cdots,{\rm d}\mathbf{z}_{\tau+H-1}(t+\frac{H-1}{H}T)]$ (9) $\displaystyle=$ $\displaystyle\tilde{f}_{\tau:\tau+H-1}(t)\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t) {\rm d}t+\tilde{g}_{\tau:\tau+H-1}(t){\rm d}\omega_{\tau:\tau+H-1}(t),$ $$

where $\tilde{f}_{\tau:\tau+H-1}(t)=[f(t),f(t+\frac{1}{H}T),\cdots,f(t+\frac{H-1}{H}T)]$ and $\tilde{g}_{\tau:\tau+H-1}(t)=[g(t),g(t+\frac{1}{H}T),\cdots,g(t+\frac{H-1}{H}T)]$. The reverse-time SDE of Eq. 9 has the following formula:

$$ $\displaystyle{\rm d}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)=$ $\displaystyle[\tilde{f}_{\tau:\tau+H-1}(t)\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)$ (10) $\displaystyle-$ $\displaystyle\tilde{g}_{\tau:\tau+H-1}(t)^{2}\nabla_{\tilde{\mathbf{z}}_{\tau: \tau+H-1}(t)}\log p_{t}\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)|\mathbf{u} _{\tau-1}(0)\big{)}]{\rm d}t$ $\displaystyle+$ $\displaystyle\tilde{g}_{\tau:\tau+H-1}(t){\rm d}\omega_{\tau:\tau+H-1}(t).$ $$

Further, the inference involving the control objective can be generalized as follows:

$$ $\displaystyle{\rm d}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)=$ $\displaystyle[\tilde{f}_{\tau:\tau+H-1}(t)\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)$ (11) $\displaystyle-$ $\displaystyle\tilde{g}_{\tau:\tau+H-1}(t)^{2}\nabla_{\tilde{\mathbf{z}}_{\tau: \tau+H-1}(t)}\log p_{t}\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)|\mathbf{u} _{\tau-1}(0)\big{)}$ $\displaystyle+$ $\displaystyle\tilde{g}_{\tau:\tau+H-1}(t)^{2}\lambda\nabla_{\tilde{\mathbf{z}} _{\tau:\tau+H-1}(t)}\mathcal{J}(\hat{\mathbf{z}}_{\tau:\tau+H-1}(0))]{\rm d}t$ $\displaystyle+$ $\displaystyle\tilde{g}_{\tau:\tau+H-1}(t){\rm d}\omega_{\tau:\tau+H-1}(t).$ $$

Here, $\hat{\mathbf{z}}_{\tau:\tau+H-1}(0)$ is an asynchronous Tweedie’s estimate (detailed in Appendix A):

$$ $\hat{\mathbf{z}}_{\tau+i}(0)=s(t+\frac{i}{H}T)^{-1}\mathbf{z}_{\tau+i}(t+\frac {i}{H}T)+s(t+\frac{i}{H}T)\sigma(t+\frac{i}{H}T)^{2}\frac{\partial\log p_{t} \big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)|\mathbf{u}_{\tau-1}(0)\big{)}}{ \partial\mathbf{z}_{\tau+i}(t+\frac{i}{H}T)}.$ $$

The distribution of $\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)$ satisfies the following proposition:

###### Proposition 1 .

Assume that the joint distribution $p\big{(}\mathbf{z}_{1}(0),\mathbf{z}_{2}(0),\cdots,\mathbf{z}_{N}(0)|\mathbf{u
}_{0}\big{)}$ has *Markov property*.
For any $t\in[0,\frac{1}{H}T]$, we have:

$$ $\displaystyle p_{t}\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)|\mathbf{u}_{ \tau-1}(0)\big{)}=\mathbb{E}_{\mathbf{z}_{\tau:\tau+H-1}(0)}[\prod_{i=0}^{H-1} p\big{(}\mathbf{z}_{\tau+i}(t+\frac{i}{H}T)|\mathbf{z}_{\tau+i}(0)\big{)}],$ (12) $$

where the expectation is about $\mathbf{z}_{\tau:\tau+H-1}(0)\sim p_{\text{\rm train}}\big{(}\mathbf{z}_{\tau:
\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0)\big{)}$, and each

$$ $p\big{(}\mathbf{z}_{\tau+i}(t+\frac{i}{H}T)|\mathbf{z}_{\tau+i}(0)\big{)}= \mathcal{N}\big{(}\mathbf{z}_{\tau+i}(t+\frac{i}{H}T);s(t+\frac{i}{H}T)\mathbf {z}_{\tau+i}(0),s(t+\frac{i}{H}T)^{2}\sigma(t+\frac{i}{H}T)^{2}\mathbf{I}\big{ )}.$ $$

Thus,
we can obtain $\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)$ for the following loss:

$$ $\mathcal{L}_{\text{asynch}}=\mathbb{E}_{\tau,t,(\mathbf{u}_{\tau-1}(0),\tilde{ \mathbf{z}}_{\tau:\tau+H-1}(t)),\bm{\epsilon}}[\|\mathbf{\bm{\epsilon}}-\bm{ \epsilon}_{\theta}\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t),\mathbf{u}_{ \tau-1},t\big{)}\|_{2}^{2}\big{]}.$ (13) $$

Here, the expectation in Eq. 13 is about $\tau\sim U(1,N-H+1)$, $t\sim U(0,\frac{1}{H}T)$, $(\mathbf{u}_{\tau-1}(0),\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t))\sim p_{\text{
train}}\big{(}\mathbf{u}_{\tau-1}(0)\big{)}p_{t}\big{(}\tilde{\mathbf{z}}_{
\tau:\tau+H-1}(t)|\mathbf{u}_{\tau-1}(0)\big{)}$, $\bm{\epsilon}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. After training, we can compute the score function by $\nabla_{\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)}\log p_{t}(\tilde{\mathbf{z}}_{
\tau:\tau+H-1}(t)|\mathbf{u}_{\tau-1}(0))\approx-\bm{\epsilon}_{\theta^{*}}
\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t),\mathbf{u}_{\tau-1}(0),t\big{)}$ with the trained parameters $\theta^{*}$.
Therefore, we can use Eq. 10, where $t$ goes from $\frac{1}{H}T$ to 0 0, to sample $\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)\sim p\big{(}\tilde{\mathbf{z}}_{\tau:
\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0),\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1
}{H}T)\big{)}$.

Figure: Algorithm 1 Closed-loop Control of CL-DiffPhyCon

###### Proposition 1 .

### 4.4 Closed-loop Control

During inference (control), our goal is not only generating control sequences and system trajectories that conform to physical dynamics but also optimizing our control objective.
Based on the two learned models, we now introduce the closed-loop control procedure $\bm{\epsilon}_{\phi}$ and $\bm{\epsilon}_{\theta}$. We first use the synchronous model $\bm{\epsilon}_{\phi}$ to produce the initial asynchronous variable $\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T)$ conditioned on the initial state $\mathbf{u}_{\text{env},0}$ by applying Eq. 5 in the horizon $[1,H]$.
Then, the online control process starts.
At each physical time step $\tau\geq 1$, we sample $\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)$ using $\mathbf{u}_{\text{env},\tau-1}$ and $\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)$ through Eq. 11.
We extract the control signal $\mathbf{w}_{\tau}(0)$ from the sampled $\mathbf{z}_{\tau}(0)$, and input the pair $(\mathbf{u}_{\text{env},\tau-1},\mathbf{w}_{\tau}(0))$ to the environment $G$, which outputs the next state $\mathbf{u}_{\text{env},\tau}$.
Then we sample a noise $\mathbf{z}_{\tau+H}(T)\sim\mathcal{N}(\mathbf{0},\sigma^{2}_{T}\mathbf{I})$, and concatenate it with the last $H-1$ components of $\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)$ to compose $\tilde{\mathbf{z}}_{\tau+1:\tau+H}(\frac{1}{H}T)$.
Now we take $\mathbf{u}_{\text{env},\tau-1}$, instead of the sampled $\mathbf{u}_{\tau-1}(0)$, and $\tilde{\mathbf{z}}_{\tau+1:\tau+H}(\frac{1}{H}T)$ to the next loop.
The whole procedure is presented in Algorithm 1 and also illustrated in Figure 2. Consequently, the joint distribution of $\{\mathbf{w}_{\tau}(0),\mathbf{u}_{\text{env},\tau}\}_{\tau=1}^{N}$ satisfies the following proposition.

###### Proposition 2 .

Using inference described in Algorithm 1, the following holds:

$$ $\displaystyle p\big{(}\mathbf{w}_{1}(0),\mathbf{u}_{\text{{\rm env}},1},\cdots ,\mathbf{w}_{N}(0),\mathbf{u}_{\text{{\rm env}},N}|\mathbf{u}_{\text{{\rm env} },0}\big{)}$ (14) $\displaystyle=\int$ $\displaystyle p_{\text{{\rm gd}}}\big{(}\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T) |\mathbf{u}_{\text{\rm env},0}\big{)}$ $\displaystyle\prod_{\tau=1}^{N}p_{\text{{\rm gd}}}\big{(}\tilde{\mathbf{z}}_{ \tau:\tau+H-1}(0)|\mathbf{u}_{\text{{\rm env}},\tau-1},\tilde{\mathbf{z}}_{ \tau:\tau+H-1}(\frac{1}{H}T)\big{)}p_{{\rm G}}\big{(}\mathbf{u}_{\text{{\rm env }},\tau}|\mathbf{u}_{\text{{\rm env}},\tau-1},\mathbf{w}_{\tau}(0)\big{)}$ $\displaystyle\mathcal{N}(\mathbf{z}_{\tau+H}(T);\mathbf{0},\sigma_{T}^{2} \mathbf{I}){\rm d}\{\mathbf{u}_{0}(0),\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T), \cdots,\mathbf{u}_{N}(0),\tilde{\mathbf{z}}_{N+1:N+H}(\frac{1}{H}T)\},$ $$

where $p_{\text{{\rm gd}}}$ denotes the transition distribution of the guided sampling (Eq. 11), and $p_{G}$ denotes the transition distribution of the environment $G$.

According to Eq. 14, for every $\tau$, the control signal $\mathbf{w}_{\tau}(0)$ is conditional on the environment states $\mathbf{u}_{\text{env},0:\tau-1}$ instead of the predicted $\mathbf{u}_{0:\tau-1}(0)$. Therefore, our method achieves closed-loop control.
It is clear that CL-DiffPhyCon is $H/h$ times faster than DiffPhyCon-$h$, the control method that conducts a full sampling process of DiffPhyCon with horizon $H$ every $h$ physical time steps, as illustrated in Figure 1. Even compared with the adaptive replaning diffusion method (Zhou et al., ), CL-DiffPhyCon is still more efficient because each latent variable is sampled only once and it does not involve extra computation such as likelihood estimation.

Additionally, although some fast sampling methods were proposed recently, such as DDIM Song et al. (), to reduce sampling cost of diffusion models, our CL-DiffPhyCon still has an independent acceleration effect beyond them. The key insight is that CL-DiffPhyCon actually does not reduce the total number of sampling steps in each physical time step, which brings the opportunity to incorporate DDIM. Hence, by introducing DDIM in both the sampling processes of the synchronous and asynchronous models, the control efficiency of CL-DiffPhyCon could be further enhanced.

###### Proposition 2 .

## 5 Experiment

In the experiments, we aim to answer three questions: (1) Can CL-DiffPhyCon outperform the classical and state-of-the-art baselines? (2) Can CL-DiffPhyCon achieve the desired acceleration of inference as we analyzed in Section 4.4, compared with diffusion control baselines, and obtain further acceleration by involving DDIM (Song et al., )? (3) Can CL-DiffPhyCon address the challenges of partial observation and high dimensional indirect control? To answer these questions, we conduct experiments on two control tasks: 1D Burgers’ equation control and 2D incompressible fluid control, both of which have important applications in science and engineering.
We provide the [code](https://anonymous.4open.science/r/CL-DiffPhyCon-C8AD) to reproduce the result for our method and baselines.

### 5.1 Baselines

The following classical, imitation learning, reinforcement learning, and diffusion control methods are selected as baselines: the classical control algorithm, widely used across various domains, Proportional-Integral-Derivative (PID) (Li et al., ); an imitation learning method Behaviour Cloning (BC) (Pomerleau, ); a recent reinforcement learning method Behavior Proximal Policy Optimization (BPPO) (Zhuang et al., ); two diffusion control methods, including RDM (Zhou et al., ), which adaptively decides when to conduct a full sampling process of control sequence, and DiffPhyCon (Wei et al., ), whose original version plans the control sequences of the whole trajectories in one sampling process. Since the trajectory length is much longer than the diffusion model horizon $H$ ($H$ is set as $H=16$ and $H=15$ in 1D and 2D tasks, respectively), we extend DiffPhyCon to its three variants DiffPhyCon-$h$ ($h\in\{1,5,H-1\}$) as our baselines, based on the predefined interval $h$ of physical time steps to conduct a full sampling process. All the diffusion control baselines use the trained synchronous diffusion model of CL-DiffPhyCon for sampling, with steps $T=900$ and $T=600$ on the 1D and 2D tasks, respectively. PID is inapplicable to the complex 2D task (Åström & Hägglund, ). RDM is reproduced following the official [code](https://github.com/rainbow979/replandiffuser). However, the default values of two thresholding hyperparameters in RDM do not well well on our tasks.
Therefore, we select a pair of values that perform best (see Appendix H for details). For other baselines, we follow the implementations in DiffPhyCon Wei et al. ().

### 5.2 1D Burgers’ Equation Control

Experiment settings. The Burgers’ equation is a widely used equation to describe a variety of physical systems. We follow the works in Hwang et al. (); Mowlavi & Nabi () and consider the 1D Burgers’ equation with the Dirichlet boundary condition and external force $w(t,x)$ as follows:

$$ $\displaystyle\begin{cases}\frac{\partial u}{\partial t}=-u\cdot\frac{\partial u }{\partial x}+\nu\frac{\partial^{2}u}{\partial x^{2}}+w(t,x)&\text{in }[0,T] \times\Omega,\\ u(t,x)=0\quad\quad\quad\quad\quad\quad\quad\quad&\text{in }[0,T]\times\partial \Omega,\\ u(0,x)=u_{0}(x)\quad\quad\quad\quad\quad\quad&\text{in }\{t=0\}\times\Omega, \end{cases}$ (15) $$

where $\nu$ is the viscosity, and $u_{0}(x)$ is the initial condition. Given a target state $u_{d}(t,x)$ defined in range $[0,T]\times\Omega$, the control objective $\mathcal{J}$ is to minimize the error between $u(t,x)$ and $u_{d}(t,x)$:

$$ $\mathcal{J}\coloneqq\int_{T}\int_{\Omega}|u(t,x)-u_{d}(t,x)|^{2}\mathrm{d}x \mathrm{d}t,$ (16) $$

subject to Eq. 15. To better meet the demands of real-world scenarios, we follow Wei et al. () and select two different experimental settings: full observation (FO) and partial observation (PO). FO means all spatial positions can be observed and controlled, while PO means only a part of them can be observed and controlled, as shown in Fig. 3.
We provide details about the dataset, experimental settings, and implementation in Appendix D.

**Table 1: Comparison results on 1D Burgers’ equation control. The average control objective $\mathcal{J}$ and inference time on a single NVIDIA A100 80GB GPU with 16 GPU cores are reported. Bold font is the best model and the runner-up is underlined.**
|  | FO | PO |  |  |
| --- | --- | --- | --- | --- |
|  | $\mathcal{J}$ $\downarrow$ | time (s) $\downarrow$ | $\mathcal{J}$ $\downarrow$ | time (s) $\downarrow$ |
| BC | 0.4708 | 0.7236 | 0.2558 | 0.7150 |
| BPPO | 0.4686 | 0.7524 | 0.2033 | 0.7342 |
| PID | 0.3250 | 0.7281 | 0.2212 | 0.7236 |
| DiffPhyCon-1 | 0.0210 | 37.1969 | 0.0196 | 39.9720 |
| DiffPhyCon-5 | 0.0252 | 8.1562 | 0.0184 | 8.4481 |
| DiffPhyCon-15 | 0.0361 | 3.7378 | 0.0192 | 3.8935 |
| RDM | 0.0296 | 6.6077 | 0.0196 | 9.8153 |
| CL-DiffPhyCon (ours) | 0.0096 | 3.4919 | 0.0090 | 3.7516 |
| CL-DiffPhyCon (DDIM, ours) | 0.0112 | 0.7585 | 0.0104 | 0.7628 |

Results. In Table 1, we present the results of our proposed CL-DiffPhyCon and baselines. Note that the reported metrics in two different settings FO/PO are not directly comparable. It can be seen that CL-DiffPhyCon delivers the best results compared to all baselines on both FO and PO settings. BC and BPPO perform poorly on this task as they rely heavily on the quality of training control sequences, which are far from optimal solutions (see Appendix D.1). The diffusion control baselines perform better, compared with BC and BPPO, because diffusion models conduct global optimization over each horizon under the guidance of control objective and excel in sampling from high dimensional space (Wei et al., ). Still, CL-DiffPhyCon decreases $\mathcal{J}$ by $54.3\%$ and $51.1\%$ compared with DiffPhyCon-1 and DiffPhyCon-5, respectively, the best baselines in the FO and PO settings. In Figure 3, we present visual comparisons of one test sample randomly selected from the test set for both FO and PO settings. They show that CL-DiffPhyCon achieves lower accumulated error regarding the target state compared with the best baselines.
The superiority of CL-DiffPhyCon over these diffusion control baselines implies that closed-loop control is necessary for good control performance, and our CL-DiffPhyCon could effectively incorporate real-time feedback from the environment into subsequent sampling of control signals.
More visualization results are presented in Appendix E.
Besides the enhanced control performance, Table 1 also reveals that CL-DiffPhyCon brings significant acceleration of the sampling process, about $H/h$ times faster than DiffPhyCon-$h$ and two times faster than RDM. Furthermore, by combining with DDIM of 30 sampling steps, CL-DiffPhyCon achieves an additional $5\times$ speedup, demonstrating the independent acceleration effect of CL-DiffPhyCon from existing fast sampling methods of diffusion models. As a result, the reduced time cost of CL-DiffPhyCon is comparable with those non-diffusion control baselines BC, BPPO, and PID, while the performance is much superior.

Figure: Figure 3: Visual comparisons between CL-DiffPhyCon (bottom row) and best baselines (middle row) on 1D Burgers’ equation control. In the figures, $J$ is calculated on the visualized samples. In the PO setting (right), white areas with oblique lines are unobservable.
Refer to caption: x3.png

**Table 2: Comparison results on 2D incompressible fluid control. The average control objective $\mathcal{J}$ and inference time on a single NVIDIA A6000 48GB GPU with 16 GPU cores are reported. Bold font is the best model and the runner-up is underlined.**
|  | Fixed map (FM) | Random map (RM) |  |  |
| --- | --- | --- | --- | --- |
|  | $\mathcal{J}$ $\downarrow$ | time (s) $\downarrow$ | $\mathcal{J}$ $\downarrow$ | time (s) $\downarrow$ |
| BC | 0.6722 | 4.61 | 0.7046 | 4.62 |
| BPPO | 0.6343 | 4.62 | 0.6524 | 4.66 |
| DiffPhyCon-1 | 0.5454 | 1677 | 0.3754 | 1692 |
| DiffPhyCon-5 | 0.5051 | 345.85 | 0.5458 | 397.95 |
| DiffPhyCon-14 | 0.5823 | 135.27 | 0.5621 | 154.28 |
| RDM | 0.4074 | 229.2 | 0.4356 | 344.4 |
| CL-DiffPhyCon (ours) | 0.3371 | 140.83 | 0.3485 | 140.54 |
| CL-DiffPhyCon (DDIM, ours) | 0.4100 | 20.75 | 0.4254 | 22.04 |

### 5.3 2D incompressible fluid control

Experiment settings.
Dynamics of 2D incompressible fluid follows the Navier-Stokes equations:

$$ $\displaystyle\begin{cases}&\frac{\partial\mathbf{v}}{\partial t}+\mathbf{v} \cdot\nabla\mathbf{v}-\nu\nabla^{2}\mathbf{v}+\nabla p=f,\\ &\nabla\cdot\mathbf{v}=0,\\ &\mathbf{v}(0,\mathbf{x})=\mathbf{v}_{0}(\mathbf{x}),\end{cases}$ (17) $$

where $f$ denotes the external force, $p$ denotes pressure, $\nu$ denotes the viscosity coefficient and $\mathbf{v}$ denotes velocity. $\mathbf{v}_{0}(\mathbf{x})$ is the initial condition. Our experimental setup is based on previous studies of Holl et al. (); Wei et al. (). Given an initial cloud of smoke in a $64\times 64$ fluid field, this task aims to minimize the volume of smoke failing to pass through the top middle exit (seven exits in total)
over $N=64$ physical time steps, by applying a sequence of 2D forces at each grid outside the outermost obstacles. This task represents a simplified scenario of many
real-life applications such as indoor air quality control (Nair et al., ).
It is very challenging for two reasons: First, the control forces can only be exerted indirectly in the peripheral regions, which necessitates the model to plan ahead across the entire trajectory to prevent the smoke from entering the wrong exits or getting stuck in corners.
Second, the dimension of control parameters in these peripheral regions is $1,792\times 2$ in each physical time step.
Therefore, this is a very challenging high-dimensional indirect control problem. In addition, we follow Zhou et al. () and add $p=0.1$ probability of random control in the execution of control signals. Random controls cause unexpected changes in the system state, which may render previously planned control sequences no longer applicable, thus necessitating re-planning and adding additional challenge to the task. Besides a fixed map (FM) setting where test samples use the same obstacles’ configuration with training trajectories, we also introduce a random map (RM) setting where the obstacles’ configuration varies in the fluid field.
For details of dataset generation, experimental settings, and implementation, please refer to Appendix F.

Figure: Figure 4: Visual comparisons between CL-DiffPhyCon and best baselines on 2D fluid control under fixed map (top) and random map (bottom) settings, respectively.
Refer to caption: x4.png

Results.
In Table 2, we show the control performance of CL-DiffPhyCon and baselines. The results indicate that CL-DiffPhyCon outperforms the baselines in both fixed map and random map settings.
This validates that our asynchronous diffusion model could effectively sample appropriate subsequent control sequences by conditioning on the changed system state in a closed-loop manner, showcasing strong generalizability.
The results also demonstrate the advantage of closed-loop diffusion control over BC and BPPO.
In Figure 4, we illustrate randomly selected test samples controlled by CL-DiffPhyCon and the best baselines. CL-DiffPhyCon exhibits a stronger capability of adapting to changed fluid and maps.
More visualization results are presented in Appendix G.
We also present the comparison of average inference time (including the time for environment execution) in Table 2.
Aligning with our observation on the 1D task, CL-DiffPhyCon achieves approximately $H/h$ times speedup compared to DiffPhyCon-$h$. Besides, CL-DiffPhyCon is much more efficient than RDM, especially in the RM setting where the larger discrepancy between training and inference necessitates more frequent replanning. Additionally, by adopting DDIM
of 75 sampling steps, the inference is further accelerated.

## 6 Conclusion

In this paper, we have proposed CL-DiffPhyCon, a novel diffusion-based method for closed-loop control of complex physical systems, based on a theoretical analysis of the target distribution. By assigning increasing noise levels within a model horizon, CL-DiffPhyCon allows for real-time interaction with the environment and efficient sampling of control signals.
Experiments on two physical control tasks demonstrate its superior control performance and efficiency.
We believe our approach is also beneficial for other tasks beyond physical systems such as robot and done control.
More discussions about limitation of our method and future work are provided in Appendix I.

## Acknowledgments

We particularly thank Chenglei Yu (Westlake University) for discussing and suggesting the solution of the diffusion models, and Zifeng Zhuang (Westlake University) for suggesting the implementation of BC and BPPO baselines. We also gratefully acknowledge the support of Westlake University Research Center for Industries of the Future and
Westlake University Center for High-performance Computing.
The content is solely the responsibility of the authors and does not necessarily represent the official views of the funding entities.

## Appendix A Basic Properties of stochastic differential equation (SDE)

In this work, we use the formulas of stochastic differential equation (SDE) to express diffusion models.
Specifically, the SDEs that we use can be written to the following unified form:

$$ ${\rm d}\mathbf{x}(t)=f(t)\mathbf{x}(t){\rm d}t+g(t){\rm d}\omega(t).$ (18) $$

The reverse-time SDE of Eq. 18 is:

$$ ${\rm d}\mathbf{x}(t)=[f(t)\mathbf{x}(t)-g(t)^{2}\nabla_{\mathbf{x}(t)}\log p_{ t}\big{(}\mathbf{x}(t)\big{)}]{\rm d}t+g(t){\rm d}\omega(t).$ (19) $$

In this section, we summarize some basic properties of Eq. 18, which can also be found in existing SDE works Song et al. (); Chung et al. (); Karras et al. ().

First, such a linear SDE has the closed solution Särkkä & Solin ()

$$ $\mathbf{x}(t)=\mathbf{x}(0)\exp(\int_{0}^{t}f(\eta){\rm d}\eta)+\int_{0}^{t} \exp(\int_{\zeta}^{t}f(\eta){\rm d}\eta)g(\zeta){\rm d}\omega(\zeta),$ (20) $$

where the second term is an Itô integral.
And the mean $\mathbf{m}(t)$ and covariance $\mathbf{V}(t)$ satisfies the following ordinary differential equations:

$$ $\displaystyle\frac{{\rm d}\mathbf{m}(t)}{{\rm d}t}$ $\displaystyle=f(t)\mathbf{m}(t),$ (21) $\displaystyle\frac{{\rm d}\mathbf{V}(t)}{{\rm d}t}$ $\displaystyle=2f(t)\mathbf{V}(t)+g(t)^{2}\mathbf{I}.$ $$

Solving these equations, we get:

$$ $\displaystyle\mathbf{m}(t)$ $\displaystyle=s(t)\mathbf{m}(0),$ (22) $\displaystyle\mathbf{V}(t)$ $\displaystyle=s(t)^{2}(\sigma(t)^{2}\mathbf{I}+\mathbf{V}(0)),$ $$

where $s(t)=\exp\big{(}\int_{0}^{t}f(\eta){\rm d}\eta\big{)}$ and $\sigma(t)=\sqrt{\int_{0}^{t}\frac{g(\eta)^{2}}{s(\eta)^{2}}{\rm d}\eta}$.

From Eq. 20 and Eq. 22, we know

$$ $\mathbf{x}(t)|\mathbf{x}(0)\sim\mathcal{N}\big{(}s(t)\mathbf{x}(0),s(t)^{2} \sigma(t)^{2}\mathbf{I}\big{)}.$ (23) $$

According to the density function of normal distribution, we have:

$$ $\nabla_{\mathbf{x}(t)}\log p\big{(}\mathbf{x}(t)|\mathbf{x}(0)\big{)}=s(t)^{-2 }\sigma(t)^{-2}(s(t)\mathbf{x}(0)-\mathbf{x}(t)).$ (24) $$

Take the expectation over $\mathbf{x}(0)$ conditional on $\mathbf{x}(t)$ on both sides:

$$ $\displaystyle s(t)^{-2}\sigma(t)^{-2}(s(t)\mathbb{E}[\mathbf{x}(0)|\mathbf{x}( t)]-\mathbf{x}(t))$ (25) $\displaystyle=$ $\displaystyle\mathbb{E}_{\mathbf{x}(0)}[\nabla_{\mathbf{x}(t)}\log p\big{(} \mathbf{x}(t)|\mathbf{x}(0)\big{)}|\mathbf{x}(t)]$ $\displaystyle=$ $\displaystyle\int\nabla_{\mathbf{x}(t)}\log p\big{(}\mathbf{x}(t)|\mathbf{x}(0 )\big{)}p\big{(}\mathbf{x}(0)|\mathbf{x}(t)\big{)}{\rm d}\mathbf{x}(0)$ $\displaystyle=$ $\displaystyle\int\nabla_{\mathbf{x}(t)}p\big{(}\mathbf{x}(t),\mathbf{x}(0)\big {)}/p\big{(}\mathbf{x}(t)\big{)}{\rm d}\mathbf{x}(0)$ $\displaystyle=$ $\displaystyle\nabla_{\mathbf{x}(t)}p\big{(}\mathbf{x}(t)\big{)}/p\big{(} \mathbf{x}(t)\big{)}$ $\displaystyle=$ $\displaystyle\nabla_{\mathbf{x}(t)}\log p\big{(}\mathbf{x}(t)\big{)}.$ $$

Rearranging the equation, we get Tweedie’s estimate:

$$ $\hat{\mathbf{x}}(0)\triangleq\mathbb{E}[\mathbf{x}(0)|\mathbf{x}(t)]=s(t)^{-1} \mathbf{x}(t)+s(t)\sigma(t)^{2}\nabla_{\mathbf{x}(t)}\log p\big{(}\mathbf{x}(t )\big{)}.$ (26) $$

Now we consider the following loss function Ho et al. ():

$$ $\displaystyle\mathbb{E}_{t,\mathbf{x}(0),\bm{\epsilon}}[\|\mathbf{\bm{\epsilon }}-\bm{\epsilon}_{\phi}\big{(}s(t)\mathbf{x}(0)+s(t)^{2}\sigma(t)^{2}\mathbf{ \bm{\epsilon}},t\big{)}\|_{2}^{2}]$ (27) $\displaystyle=$ $\displaystyle\mathbb{E}_{t,\mathbf{x}(0),\mathbf{x}(t)}[\|\frac{\mathbf{x}(t)- \mathbf{x}(0)s(t)}{s(t)^{2}\sigma(t)^{2}}-\bm{\epsilon}_{\phi}\big{(}\mathbf{x }(t),t\big{)}\|_{2}^{2}].$ $$

The best prediction $\bm{\epsilon}_{\phi^{*}}(\cdot,t)$ in this loss is the following conditional expectation:

$$ $\displaystyle\bm{\epsilon}_{\phi^{*}}\big{(}\mathbf{x}(t),t\big{)}=$ $\displaystyle\mathbb{E}[\frac{\mathbf{x}(t)-\mathbf{x}(0)s(t)}{s(t)^{2}\sigma( t)^{2}}-\bm{\epsilon}_{\phi}\big{(}\mathbf{x}(t),t\big{)}|\mathbf{x}(t)]$ (28) $\displaystyle=$ $\displaystyle\frac{\mathbf{x}(t)-\mathbb{E}[\mathbf{x}(0)|\mathbf{x}(t)]s(t)}{ s(t)^{2}\sigma(t)^{2}}.$ $$

Plugging Eq. 26 in it, we get the score function:

$$ $\displaystyle\nabla_{\mathbf{x}(t)}\log p\big{(}\mathbf{x}(t)\big{)}=-\bm{ \epsilon}_{\phi^{*}}\big{(}\mathbf{x}(t),t\big{)}.$ (29) $$

## Appendix B Variance preserving (VP) SDE

In the practical training and inference, we use Variance Preserving (VP) SDE Song et al. (), which is actually the continous version of DDPM Ho et al. ().

VP SDE specializes Eq. 18 to

$$ ${\rm d}\mathbf{x}(t)=-\frac{1}{2}\beta(t)\mathbf{x}(t){\rm d}t+\sqrt{\beta(t)} {\rm d}\omega(t),$ (30) $$

where $\beta(t)>0$ for $t\in[0,T]$.
And its reverse-time SDE is:

$$ ${\rm d}\mathbf{x}(t)=[-\frac{1}{2}\beta(t)\mathbf{x}(t)-\beta(t)\nabla_{ \mathbf{x}(t)}\log p_{t}\big{(}\mathbf{x}(t)\big{)}]{\rm d}t+\sqrt{\beta(t)}{ \rm d}\omega(t),$ (31) $$

In algorithm implementation, we use $K$ time steps, $t_{i}=\frac{i}{K}T$ for $i=0,\cdots,K-1$, to discretize Eq. 30.
Using $\Delta t=\frac{1}{K}T$, we have:

$$ $\displaystyle\mathbf{x}(t_{i}+\Delta t)$ $\displaystyle\approx\mathbf{x}(t_{i})-\frac{1}{2}\beta(t_{i})\Delta t\mathbf{x }(t_{i})+\sqrt{\beta(t_{i})\Delta t}\xi,$ (32) $\displaystyle\approx\sqrt{1-\beta(t_{i})\Delta t}\mathbf{x}(t_{i})+\sqrt{\beta (t_{i})\Delta t}\xi,$ $$

where $\xi\sim\mathcal{N}(\mathbf{0},\mathbf{I})$.

Let $\beta_{i}=\beta(t_{i})\Delta t$, $\mathbf{x}_{i}=\mathbf{x}(t_{i})$.
We get the discrete forward iterative equation, which is exact DDPM:

$$ $\displaystyle\mathbf{x}_{i+1}=\sqrt{1-\beta_{i}}\mathbf{x}_{i}+\sqrt{\beta_{i} }\xi.$ (33) $$

Similarly, Eq. 31 can be discretize to

$$ $\displaystyle\mathbf{x}(t_{i}-\Delta t)$ $\displaystyle\approx\mathbf{x}(t_{i})+\frac{1}{2}\beta(t_{i})\Delta t\mathbf{x }(t_{i})+\beta(t_{i})\Delta t\nabla_{\mathbf{x}(t_{i})}\log p_{t_{i}}\big{(} \mathbf{x}(t_{i})\big{)}+\sqrt{\beta(t_{i})\Delta t}\xi$ (34) $\displaystyle\approx\frac{1}{\sqrt{1-\beta(t_{i})\Delta t}}\mathbf{x}(t_{i})+ \beta(t_{i})\Delta t\nabla_{\mathbf{x}(t_{i})}\log p_{t_{i}}\big{(}\mathbf{x}( t_{i})\big{)}+\sqrt{\beta(t_{i})\Delta t}\xi.$ $$

So, we get the discrete backward iterative equation:

$$ $\displaystyle\mathbf{x}_{i-1}=\frac{1}{\sqrt{1-\beta_{i}}}\mathbf{x}_{i}+\beta _{i}\nabla_{\mathbf{x}_{i}}\log p\big{(}\mathbf{x}_{i}\big{)}+\sqrt{\beta_{i}}\xi.$ (35) $$

Additionally, we have:

$$ $\displaystyle s(t_{i})^{2}$ $\displaystyle=\exp\big{(}\int_{t_{0}}^{t_{i}+\Delta t}-\beta(\eta){\rm d}\eta \big{)}$ (36) $\displaystyle=\exp\big{(}\sum_{j=0}^{i}-\beta(t_{j})\Delta t\big{)}$ $\displaystyle=\prod_{j=0}^{i}\exp\big{(}-\beta(t_{j})\Delta t\big{)}$ $\displaystyle\approx\prod_{j=0}^{i}(1-\beta(t_{j})\Delta t)$ $\displaystyle=\prod_{j=0}^{i}(1-\beta_{j}).$ $$

And, we also have:

$$ $\displaystyle s(t_{i})^{2}\sigma(t_{i})^{2}$ $\displaystyle=\exp\big{(}\int_{t_{0}}^{t_{i}+\Delta t}-\beta(\eta){\rm d}\eta \big{)}\int_{t_{0}}^{t_{i}+\Delta t}\beta(\eta)\exp\big{(}\int_{t_{0}}^{\eta} \beta(\zeta){\rm d}\zeta\big{)}{\rm d}\eta$ (37) $\displaystyle=\int_{t_{0}}^{t_{i}+\Delta t}\beta(\eta)\exp\big{(}\int_{\eta}^{ t_{i}+\Delta t}-\beta(\zeta){\rm d}\zeta\big{)}{\rm d}\eta$ $\displaystyle\approx 1-\exp\big{(}\int_{t_{0}}^{t_{i}+\Delta t}-\beta(\eta) \exp\big{(}\int_{\eta}^{t_{i}+\Delta t}-\beta(\zeta){\rm d}\zeta\big{)}{\rm d} \eta\big{)}$ $\displaystyle\approx 1-\exp\big{(}\int_{t_{0}}^{t_{i}+\Delta t}-\beta(\eta){ \rm d}\eta\big{)}$ $\displaystyle=1-s(t_{i})^{2}$ $\displaystyle\approx 1-\prod_{j=0}^{i}(1-\beta_{j}).$ $$

We define $\alpha_{i}=1-\beta_{i}$ and $\bar{\alpha}_{i}=\prod_{j=0}^{i}\alpha_{i}$.
Then, we get $s(t_{i})^{2}\approx\bar{\alpha}_{i}$, $s(t_{i})^{2}\sigma(t_{i})^{2}\approx 1-\bar{\alpha}_{i}$, and $\mathbf{x}_{i}|\mathbf{x}_{0}\sim\mathcal{N}\big{(}\sqrt{\bar{\alpha}_{i}}
\mathbf{x}_{0},(1-\bar{\alpha}_{i})\mathbf{I}\big{)}$.
Therefore, Tweedie’s estimate becomes

$$ $\hat{\mathbf{x}}(0)=\frac{1}{\sqrt{\bar{\alpha}_{i}}}(x_{i}+(1-\bar{\alpha}_{i })\nabla_{\mathbf{x}_{i}}\log p(\mathbf{x}_{i})).$ (38) $$

According to Eq. 22, the covariance of $\mathbf{x}_{i}$ is $(1-\bar{\alpha}_{i})\mathbf{I}+\bar{\alpha}_{i}\mathbf{V}_{0}$.
Assuming the clear data has covariance $\mathbf{I}$, we can start from $\mathbf{x}_{K}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ during inference.

## Appendix C Derivations

###### Proof of Theorem 1 .

We first conduct the following decomposition:

$$ $\displaystyle p\big{(}\mathbf{z}_{1:N}(0),\tilde{\mathbf{z}}_{1:H}(\frac{1}{H} T),\cdots,\tilde{\mathbf{z}}_{N+1:N+H}(\frac{1}{H}T)|\mathbf{u}_{0}\big{)}$ (39) $\displaystyle=$ $\displaystyle p\big{(}\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T),\mathbf{z}_{1}(0) ,\tilde{\mathbf{z}}_{2:H+1}(\frac{1}{H}T),\cdots,\mathbf{z}_{N}(0),\tilde{ \mathbf{z}}_{N+1:N+H}(\frac{1}{H}T)|\mathbf{u}_{0}\big{)}$ $\displaystyle=$ $\displaystyle p\big{(}\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T)|\mathbf{u}_{0} \big{)}$ $\displaystyle\prod_{\tau=1}^{N}p\big{(}\mathbf{z}_{\tau}(0),\tilde{\mathbf{z}} _{\tau+1:\tau+H}(\frac{1}{H}T)|\mathbf{u}_{0},\tilde{\mathbf{z}}_{1:H}(\frac{1 }{H}T),\mathbf{z}_{1}(0),\tilde{\mathbf{z}}_{2:H+1}(\frac{1}{H}T),\cdots, \mathbf{z}_{\tau-1}(0),\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}$ $\displaystyle=$ $\displaystyle p\big{(}\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T)|\mathbf{u}_{0} \big{)}$ $\displaystyle\prod_{\tau=1}^{N}p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0), \mathbf{z}_{\tau+H}(T)|\mathbf{u}_{0},\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T), \mathbf{z}_{1}(0),\tilde{\mathbf{z}}_{2:H+1}(\frac{1}{H}T),\cdots,\mathbf{z}_{ \tau-1}(0),\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}$ $\displaystyle=$ $\displaystyle p\big{(}\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T)|\mathbf{u}_{0} \big{)}\prod_{\tau=1}^{N}\mathcal{N}(\mathbf{z}_{\tau+H}(T)|\mathbf{0},\sigma_ {T}^{2}\mathbf{I})$ $\displaystyle\prod_{\tau=1}^{N}p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)| \mathbf{u}_{0},\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T),\mathbf{z}_{1}(0),\tilde {\mathbf{z}}_{2:H+1}(\frac{1}{H}T),\cdots,\mathbf{z}_{\tau-1}(0),\tilde{ \mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}.$ $$

The third equation holds since $[\mathbf{z}_{\tau}(0),\tilde{\mathbf{z}}_{\tau+1:\tau+H}(\frac{1}{H}T)]$ and $[\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0),\mathbf{z}_{\tau+H}(T)]$ are two different arrangements of a same vector $[\mathbf{z}_{\tau}(0),\mathbf{z}_{\tau+1}(\frac{1}{H}T),\cdots,\mathbf{z}_{
\tau+H-1}(\frac{H-1}{H}T),\mathbf{z}_{\tau+H}(T)]$.
The last equation holds since $\mathbf{z}_{\tau}(T)$ is independently normally distributed with density $\mathcal{N}(\mathbf{z}_{\tau}(T)|\mathbf{0},\sigma_{T}^{2}\mathbf{I})$ for any $\tau$. Then, we analyze the conditional probability $p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{0},\tilde{\mathbf{z}
}_{1:H}(\frac{1}{H}T),\mathbf{z}_{1}(0),\tilde{\mathbf{z}}_{2:H+1}(\frac{1}{H}
T),\cdots,\mathbf{z}_{\tau-1}(0),\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H
}T)\big{)}$ in detail.

Due to Markov property of physical systems, when $\mathbf{u}_{\tau-1}(0)$ appears, other variables with temporal subscripts less than $\tau$ can not affect those with subscripts greater than or equal to $\tau$.
Therefore, we can remove those variables with subscripts less than $\tau$, except for $\mathbf{u}_{\tau-1}(0)$. Thus, we have:

$$ $\displaystyle p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{0}, \tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T),\mathbf{z}_{1}(0),\tilde{\mathbf{z}}_{2 :H+1}(\frac{1}{H}T),\cdots,\mathbf{z}_{\tau-1}(0),\tilde{\mathbf{z}}_{\tau: \tau+H-1}(\frac{1}{H}T)\big{)}$ (40) $\displaystyle=$ $\displaystyle p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1 }(0),\mathbf{c},\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)},$ $$

where $\mathbf{c}$ denotes the subset of latent variables $[\mathbf{z}_{\tau}(T),\tilde{\mathbf{z}}_{\tau:\tau+1}(\frac{H-1}{H}T),\cdots,
\tilde{\mathbf{z}}_{\tau:\tau+H-2}(\frac{2}{H}T)]$ extracted from $[\tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T),\cdots,\mathbf{z}_{\tau-1}(0),\tilde{
\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)]$. Further, by adding new variables $\mathbf{z}_{\tau:\tau+H-1}(0)$ to this probability, we have

$$ $\displaystyle p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1 }(0),\mathbf{c},\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}$ (41) $\displaystyle=$ $\displaystyle\int p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{ \tau-1}(0),\mathbf{c},\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T),\mathbf {z}_{\tau:\tau+H-1}(0)\big{)}$ $\displaystyle p\big{(}\mathbf{z}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0), \mathbf{c},\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}{\rm d} \mathbf{z}_{\tau:\tau+H-1}(0).$ $$

The first distribution inside the integral in Eq. 41 can be simplified as follows:

$$ $\displaystyle p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1 }(0),\mathbf{c},\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T),\mathbf{z}_{ \tau:\tau+H-1}(0)\big{)}$ (42) $\displaystyle=$ $\displaystyle\frac{p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0),\mathbf{c}, \tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T),\mathbf{z}_{\tau:\tau+H-1}(0) |\mathbf{u}_{\tau-1}(0)\big{)}}{p\big{(}\mathbf{c},\tilde{\mathbf{z}}_{\tau: \tau+H-1}(\frac{1}{H}T),\mathbf{z}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0) \big{)}}$ $\displaystyle=$ $\displaystyle\frac{p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0),\tilde{ \mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T),\mathbf{z}_{\tau:\tau+H-1}(0)| \mathbf{u}_{\tau-1}(0)\big{)}p\big{(}\mathbf{c}|\tilde{\mathbf{z}}_{\tau:\tau+ H-1}(\frac{1}{H}T)\big{)}}{p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1} {H}T),\mathbf{z}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0)\big{)}p\big{(} \mathbf{c}|\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}}$ $\displaystyle=$ $\displaystyle\frac{p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0),\tilde{ \mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T),\mathbf{z}_{\tau:\tau+H-1}(0)| \mathbf{u}_{\tau-1}(0)\big{)}}{p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}( \frac{1}{H}T),\mathbf{z}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0)\big{)}}$ $\displaystyle=$ $\displaystyle p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1 }(0),\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T),\mathbf{z}_{\tau:\tau+H- 1}(0)\big{)}.$ $$

In the second equation, we use the Markov property of the diffusion process of $\mathbf{z}_{\tau}(t)$ over $t$, and skip a step including the following two equations

$$ $\displaystyle p\big{(}\mathbf{c}|\mathbf{u}_{\tau-1}(0),\mathbf{z}_{\tau:\tau+ H-1}(0),\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0),\tilde{\mathbf{z}}_{\tau:\tau+H- 1}(\frac{1}{H}T)\big{)}$ $\displaystyle=p\big{(}\mathbf{c}|\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H }T)\big{)}$ $\displaystyle p\big{(}\mathbf{c}|\mathbf{u}_{\tau-1}(0),\mathbf{z}_{\tau:\tau+ H-1}(0),\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}$ $\displaystyle=p\big{(}\mathbf{c}|\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H }T)\big{)}$ $$

in the numerator and denominator respectively of the right side of the second equation. Similarly, we show that we can also simplify the other distribution inside the integral in Eq. 41.

$$ $\displaystyle p\big{(}\mathbf{z}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0), \mathbf{c},\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}$ (43) $\displaystyle=$ $\displaystyle\frac{p\big{(}\mathbf{z}_{\tau:\tau+H-1}(0),\mathbf{c},\tilde{ \mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)|\mathbf{u}_{\tau-1}(0)\big{)}}{p\big {(}\mathbf{c},\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)|\mathbf{u}_{ \tau-1}(0)\big{)}}$ $\displaystyle=$ $\displaystyle\frac{p\big{(}\mathbf{z}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0 )\big{)}p\big{(}\mathbf{c},\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)| \mathbf{u}_{\tau-1}(0),\mathbf{z}_{\tau:\tau+H-1}(0)\big{)}}{p\big{(}\mathbf{c },\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)|\mathbf{u}_{\tau-1}(0)\big{ )}}$ $\displaystyle=$ $\displaystyle\frac{p\big{(}\mathbf{z}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0 )\big{)}p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)|\mathbf{u}_{ \tau-1}(0),\mathbf{z}_{\tau:\tau+H-1}(0)\big{)}p\big{(}\mathbf{c}|\tilde{ \mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}}{p\big{(}\tilde{\mathbf{z}}_{ \tau:\tau+H-1}(\frac{1}{H}T)|\mathbf{u}_{\tau-1}(0)\big{)}p\big{(}\mathbf{c}| \tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}}$ $\displaystyle=$ $\displaystyle\frac{p\big{(}\mathbf{z}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0 )\big{)}p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)|\mathbf{u}_{ \tau-1}(0),\mathbf{z}_{\tau:\tau+H-1}(0)\big{)}}{p\big{(}\tilde{\mathbf{z}}_{ \tau:\tau+H-1}(\frac{1}{H}T)|\mathbf{u}_{\tau-1}(0)\big{)}}$ $\displaystyle=$ $\displaystyle p\big{(}\mathbf{z}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0), \tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}.$ $$

Combining Eq. 40 to Eq. 43, we get the following result:

$$ $\displaystyle p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{0}, \tilde{\mathbf{z}}_{1:H}(\frac{1}{H}T),\mathbf{z}_{1}(0),\tilde{\mathbf{z}}_{2 :H+1}(\frac{1}{H}T),\cdots,\mathbf{z}_{\tau-1}(0),\tilde{\mathbf{z}}_{\tau: \tau+H-1}(\frac{1}{H}T)\big{)}$ (44) $\displaystyle=$ $\displaystyle\int p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{ \tau-1}(0),\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T),\mathbf{z}_{\tau: \tau+H-1}(0)\big{)}$ $\displaystyle p\big{(}\mathbf{z}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0), \tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}{\rm d}\mathbf{z}_{\tau :\tau+H-1}(0)$ $\displaystyle=$ $\displaystyle\int p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0),\mathbf{z}_{ \tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0),\tilde{\mathbf{z}}_{\tau:\tau+H-1}( \frac{1}{H}T)\big{)}{\rm d}\mathbf{z}_{\tau:\tau+H-1}(0)$ $\displaystyle=$ $\displaystyle p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1 }(0),\tilde{\mathbf{z}}_{\tau:\tau+H-1}(\frac{1}{H}T)\big{)}.$ $$

Plugging this to Eq. 39, we prove the conclusion.
∎

###### Proof of Proposition 1 .

Adding new variables $\mathbf{z}_{\tau:\tau+H-1}(0)$, we have:

$$ $\displaystyle p_{t}\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)|\mathbf{u}_{ \tau-1}(0)\big{)}$ (45) $\displaystyle=$ $\displaystyle\int p\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)|\mathbf{u}_{ \tau-1}(0),\mathbf{z}_{\tau:\tau+H-1}(0)\big{)}p\big{(}\mathbf{z}_{\tau:\tau+H -1}(0)|\mathbf{u}_{\tau-1}(0)\big{)}{\rm d}\mathbf{z}_{\tau:\tau+H-1}(0).$ $$

Notice that the distributions of $\{\mathbf{z}_{\tau+i}(t+\frac{i}{H}T)\}_{i=0}^{H-1}$ are conditional independent.
Therefore, we can reformulate the above equation to:

$$ $\displaystyle p_{t}\big{(}\tilde{\mathbf{z}}_{\tau:\tau+H-1}(t)|\mathbf{u}_{ \tau-1}(0)\big{)}$ (46) $\displaystyle=$ $\displaystyle\int\prod_{i=0}^{H-1}p\big{(}\mathbf{z}_{\tau+i}(t+\frac{i}{H}T)| \mathbf{u}_{\tau-1}(0),\mathbf{z}_{\tau:\tau+H-1}(0)\big{)}p\big{(}\mathbf{z}_ {\tau:\tau+H-1}(0)|\mathbf{u}_{\tau-1}(0)\big{)}{\rm d}\mathbf{z}_{\tau:\tau+H -1}(0)$ $\displaystyle=$ $\displaystyle\int\prod_{i=0}^{H-1}p\big{(}\mathbf{z}_{\tau+i}(t+\frac{i}{H}T)| \mathbf{z}_{\tau+i}(0)\big{)}p\big{(}\mathbf{z}_{\tau:\tau+H-1}(0)|\mathbf{u}_ {\tau-1}(0)\big{)}{\rm d}\mathbf{z}_{\tau:\tau+H-1}(0)$ $\displaystyle=$ $\displaystyle\mathbb{E}_{\mathbf{z}_{\tau:\tau+H-1}(0)}[\prod_{i=0}^{H-1}p\big {(}\mathbf{z}_{\tau+i}(t+\frac{i}{H}T)|\mathbf{z}_{\tau+i}(0)\big{)}].$ $$

Here, the second equation is because $\mathbf{z}_{\tau+i}(t+\frac{i}{H}T)$ only depends on $\mathbf{z}_{\tau+i}(0)$ with the distribution Eq. 23 described in Section A.
And we have proved the desired conclusion.
∎

###### Proof of Proposition 2 .

This proposition is just plugging the $p_{{\rm G}}\big{(}\mathbf{u}_{\text{{\rm env}},\tau}|\mathbf{u}_{\text{{\rm env
}},\tau-1},\mathbf{w}_{\tau}(0)\big{)}$ to the joint distribution in Proposition 1.
And we omit this trivial proof.
∎

###### Proof of Theorem 1 .

###### Proof of Proposition 1 .

###### Proof of Proposition 2 .

## Appendix D Details of 1D Burgers’ equation control

### D.1 Dataset

We follow instructions in Wei et al. () to generate a 1D Burgers’ equation dataset. Specifically, we use the finite difference method (FDM) to generate trajectories in a domain of space range $x\in[0,1]$ and time range $t\in[0,1]$, with random initial states and control sequences following certain distributions. The space is discretized into 128 grids and time into 10000 steps. We generated 90000 trajectories for the training set and 50 for the testing set.
For each training sample, its target state $u_{d}(t,x)$ is randomly selected from other training samples, which means that almost all the samples in the training set are unsuccessful as they could hardly achieve the target states under random controls. Thus it is challenging to generate control sequences with performance beyond the training dataset.

### D.2 Experimental Settings

Similar to (Wei et al., ), we design two settings of 1D Burgers’ equation control in Section 5.2:

Full observation: This is a common scenario where all states $u(t,x),x\in[0,1]$ for $t\in[0,1]$ of the system can be observed.

Partial observation: In this setting, we hide some parts of $\mathbf{u}$ and measure the $\mathcal{J}$ of model control as the white area with oblique lines shown in Figure 3.
Specifically, $\mathbf{u}(t,x),~{}x\in[\frac{1}{4},\frac{3}{4}]$ is set to zero in the dataset during training and $\mathbf{u}_{0}(x),~{}x\in[\frac{1}{4},\frac{3}{4}]$ is also set to zero during testing. Only $\Omega=[0,\frac{1}{4}]\cup[\frac{3}{4},1]$ is observed, controlled and evaluated. This setting is particularly challenging because of the uncertainty introduced by the unobserved states.

**Table 3: Hyperparameters of models and training for 1D Burgers’ equation control.**
| Hyperparameter name | Full observation | Partial observation |
| --- | --- | --- |
| U-Net $\bm{\epsilon}_{\phi}(\mathbf{w})$ |  |  |
| Model horizon $H$ | 16 | 16 |
| Initial dimension | 64 | 64 |
| Downsampling/Upsampling layers | 4 | 4 |
| Convolution kernel size | 3 | 3 |
| Dimension multiplier | $[1,2,4,8]$ | $[1,2,4,8]$ |
| Attention hidden dimension | 32 | 32 |
| Attention heads | 4 | 4 |
| U-Net $\bm{\epsilon}_{\theta}(\mathbf{u},\mathbf{w})$ |  |  |
| Model horizon $H$ | 16 | 16 |
| Initial dimension | 64 | 64 |
| Downsampling/Upsampling layers | 4 | 4 |
| Convolution kernel size | 3 | 3 |
| Dimension multiplier | $[1,2,4,8]$ | $[1,2,4,8]$ |
| Attention hidden dimension | 32 | 32 |
| Attention heads | 4 | 4 |
| Training |  |  |
| Training batch size | 16 | 16 |
| Optimizer | Adam | Adam |
| Learning rate | 1e-4 | 1e-4 |
| Training steps | 190000 | 190000 |
| Learning rate scheduler | cosine annealing | cosine annealing |
| Inference |  |  |
| Synchronously sampling steps | 900 | 900 |
| Each asynchronously sampling step | 60 | 60 |

### D.3 Implementation

We use U-Net (Ronneberger et al., ) as architectures for the models $\bm{\epsilon}_{\phi}$ and $\bm{\epsilon}_{\theta}$.
Two models are separately trained using the same training dataset. Note that in the partial observation settings, the unobserved data is invisible to the model during both training and testing as introduced in Appendix D.2. We simply pad zero in the corresponding locations of the model input and conditions, and also exclude these locations in the training loss. Therefore, the model only learns the correlation between the observed states and control sequences. We use the MSE loss to train the models. Both models have $T=900$ diffusion steps. The DDIM (Song et al., ) sampling we use only has 30 diffusion steps with the hyperparameter $\eta=1$. Hyperparameters of models and training are listed in Table 3.

## Appendix E 1D visualization results

We present the visualization results of our method and baselines under two settings: FO and PO in Figure 5, 6, 7, and 8, respectively. Under each setting, we present the results of four randomly selected samples from the test set.

Figure: Figure 5: Two visualizations results of 1D Burgers’ equation control under the FO (full observation) setting. The first line is the target $u_{d}(t,x)$ and the error $u(t,x)-u_{d}(t,x)$ measures the gap between the state under control and target. The horizontal axis is the time coordinate and the vertical axis is the state.
Refer to caption: x5.png

Figure: Figure 6: Two visualizations results of 1D Burgers’ equation control under the FO setting. The first line is the target $u_{d}(t,x)$ and the error $u(t,x)-u_{d}(t,x)$ measures the gap between the state under control and target. The horizontal axis is the time coordinate and the vertical axis is the state.
Refer to caption: x7.png

Figure: Figure 7: Two visualizations results of 1D Burgers’ equation control under the PO (partial observation) setting. The first line is the target $u_{d}(t,x)$ and the error $u(t,x)-u_{d}(t,x)$ measures the gap between the state under control and target. The horizontal axis is the time coordinate and the vertical axis is the state. The white area with oblique lines in the middle represents an unobservable area.
Refer to caption: x9.png

Figure: Figure 8: Two visualizations results of 1D Burgers’ equation control under the PO setting. The first line is the target $u_{d}(t,x)$ and the error $u(t,x)-u_{d}(t,x)$ measures the gap between the state under control and target. The horizontal axis is the time coordinate and the vertical axis is the state. The white area with oblique lines in the middle represents an unobservable area.
Refer to caption: x11.png

From these visualization results, it can be observed that under the control of our proposed CL-DiffPhyCon, the state $\mathbf{u}(t,x)$ is more close to the target state $\mathbf{u}_{d}(t,x)$ given different initial states. Furthermore, this observation is consistent under both settings: FO and PO, implying that our CL-DiffPhyCon is effective in addressing the partial observation challenge.

**Table 4: Hyperparameters of 2D experiments.**
| Hyperparameter Name | Value |
| --- | --- |
| Number of attention heads | 4 |
| Kernel size of conv3d | (3, 3, 3) |
| Padding of conv3d | (1,1,1) |
| Stride of conv3d | (1,1,1) |
| Kernel size of downsampling | (1, 4, 4) |
| Padding of downsampling | (1, 2, 2) |
| Stride of downsampling | (0, 1, 1) |
| Kernel size of upsampling | (1, 4, 4) |
| Padding of upsampling | (1, 2, 2) |
| Stride of upsampling | (0, 1, 1) |

## Appendix F Details of 2D incompressible fluid control

### F.1 Dataset

We follow instructions in (Wei et al., ) to generate a 2D incompressible fluid control dataset, where we use the Phiflow solver Holl et al. () to simulate the fluid dynamics. The resolution of the 2D flow field is set to be 64 $\times$ 64. The flow field is set to be boundless. The locations of obstacles are shown in the top two rows in Figure 4, which are shared for all the training trajectories. Each trajectory includes $N=64$ physical time steps of the following features: horizontal and vertical velocities, smoke density, and horizontal and vertical control forces. In total, we generated 40,000 training trajectories.

### F.2 Experimental setting

Following RDM (Zhou et al., ), we add a small level of control noise to make the task more challenging. Specifically, for a trajectory of length $N=64$, the control signal in each physical time step is executed with probability $p=0.1$ as a random control, where the horizontal and vertical components following the uniform distributions bounded by the minimal and maximal values of horizontal and vertical control signals in the training dataset, respectively.

In addition, we design two settings of fixed map (FM) and random map (RM) to test the generalization capability of each method:

Fixed map (FM): In this setting, all 50 test samples use the same obstacle configuration with training trajectories. However, the initial locations of test samples are different.

Random map (RM): In this setting, the 50 test samples use random obstacles’ configuration in the fluid field. Specifically, the movement ranges for the five internal obstacles, whose default locations are shown in the top two rows in Figure 4, are as follows:

- •
The two upper obstacles can move downward, left, or right by no more than 3 grid spaces.
- •
The two middle obstacles can move upward, downward, left, or right by no more than 3 grid spaces.
- •
The one lower obstacle can move upward, left, or right by no more than 3 grid spaces.

Note that also the obstacles’ configuration varies for different test samples, for each test sample, its obstacles’ configuration remains unchanged during the control process.

### F.3 Implementation

We use 3D U-net Ho et al. () as architectures for the models $\bm{\epsilon}_{\phi}$ and $\bm{\epsilon}_{\theta}$. Two models are separately trained using the same training dataset. We use the MSE loss to train the models. Both models have $T=600$ diffusion steps. The DDIM (Song et al., ) sampling we use has 75 diffusion steps with the hyperparameter $\eta=0.3$. Hyperparameters of models and training are listed in Table 4.

## Appendix G 2D visualization results

More visualization results of our method and compared baselines are presented.
For both fixed map (Figure 9 and Figure 10) and random map (Figure 11 and Figure 12) settings, we present two randomly selected test samples. Our method outperforms baselines on average. Note that the fluid dynamics is very sensitive to the generated control signals, thus different methods may perform dramatically differently on the same test sample. However, we find that each method performs stably on each test sample through multiple evaluations.

Figure: Figure 9: Visualizations results of 2D fluid control by our CL-DiffPhyCon method and baselines for a same test sample. This is an example of the *fixed map* setting. Each row shows six frames of smoke density. The smaller the value of $J$, the better the performance.
Refer to caption: x13.png

Figure: Figure 10: Visualizations results of 2D fluid control by our CL-DiffPhyCon method and baselines for a same test sample. This is an example of the *fixed map* setting. Each row shows six frames of smoke density. The smaller the value of $J$, the better the performance.
Refer to caption: x14.png

Figure: Figure 11: Visualizations results of 2D fluid control by our CL-DiffPhyCon method and baselines for a same test sample. This is an example of the *random map* setting. Each row shows six frames of smoke density. The smaller the value of $J$, the better the performance.
Refer to caption: x15.png

Figure: Figure 12: Visualizations results of 2D fluid control by our CL-DiffPhyCon method and baselines for a same test sample. This is an example of the *random map* setting. Each row shows six frames of smoke density. The smaller the value of $J$, the better the performance.
Refer to caption: x16.png

## Appendix H Baselines

For RDM Zhou et al. (), it involves two hyperparameters to determine when to update the previously planned control sequences or replan from scratch.
Their default values do not work well on our 1D and 2D tasks. On the 1D task, they are too small, which results in very expensive replanning from scratch almost every physical time, although its performance is still unsatisfactory. On the 2D task, they are too large, leading to few replannings and inferior results.
Therefore, we run extensive evaluations to search for their best values.
On the 1D task, we find that the best pair of values for the two hyperparameters is (0.06, 0.1) in the full observation (FO) setting and (0.15, 0.3) in the partial observation (PO) setting.
On the 2D task, we find that the best pair of values is (0.0005, 0.002).

## Appendix I Limitation and Future Work

There are still several limitations of CL-DiffPhyCon that provide exciting opportunities for future work. Firstly, CL-DiffPhyCon is currently trained offline, without interacting with the environment. Integrating the environment into the training process could allow the model to dynamically adapt to the environment and uncover new strategies and solutions.
Secondly, CL-DiffPhyCon learns the diffusion constraints directly from observed physical data.
Another clearly beneficial approach is to incorporate other physical priors and constraints, such as the formula of PDEs, into the training of diffusion models to enable better and faster learning of physical constraints. Finally, the effectiveness of CL-DiffPhyCon in more domains is worth exploring. Although our method is designed for the control of complex physical systems, it holds promise for broader applications, such as robot control and drone control.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: Detailed comments are summarized below:

1. How restrictive is the consideration of noise-free control? There are many noise issues arising from complex engineering systems and applications. Can the proposed method be extended to account for plant, control, and measurement noise? 

2. Another important aspect associated with complex and practical physical systems is physical constraints. For example, physical control actions are implemented via actuators which often have upper and lower limits. That is one major reason model predictive control is useful. Can the proposed method handle physical constraints associated with control actions? 

3.  From Figure 2 and the limitation mentioned in the appendix, the algorithm is trained offline. In that case, it is not clear how the control is designed and employed in a real environment, as controllers are typically designed and used in practice in an interactive manner. In Figure 2, the ending time instance is T, can it be changed to T+1, as a bigger loop, to make it useful in an online/interactive manner? 

4. In both examples, the controls are defined as a spatiotemporal function in the PDE model, which is too ideal to be implemented in a realistic system. For PDE system control, typically, the control is implemented as a temporal function. The corresponding actuator is located in the boundary or a small part of the entire spatial domain. It would be interesting to see the results in these more realistic settings by using the proposed method. 

5. In the Burgers’ equation example, there is no difference between the full observation and partial observation cases, because the latter case only observes and controls two parts of the entire domain. In the case considered in the manuscript, one can simply treat each observed sub-domain as a total entire domain and apply the algorithm on the PDE with new spatial coordinates/ranges/boundary conditions. Partial observation typically refers to the case using partial observation to control the entire domain. The control performance results of the unobserved state should be added to fully demonstrate the effectiveness of the proposed algorithm. Similarly. practical applications often have sensors installed in the boundary of PDE systems. It would be good to also add results by using boundary observation data. 

6. The limitations of the proposed method should be moved to the conclusion section for better readability.
soundness: 3
strengths: One advantage is to improve the computation efficiency, which is achieved by an asynchronous denoising framework. This aspect is validated using numerical simulations. Overall, this manuscript is well-written and organized.
summary: This manuscript addresses a closed-loop diffusion control problem of complex physical systems. The proposed CL-DiffPhyCon (Wei et al., 2024) extends the work of DiffPhyCon (Wei et al., 2024).
weaknesses: There are some technical limitations, including noise-free inputs, input constraints, and numerical validations by using full in-domain control and observation.

---

confidence: 4
contribution: 4
initial_rating: 10
presentation: 4
questions: I am not quite sure why the proposed method achieves a better control objective than DiffPhyCon-1 - isn't the latter replanning at each control step by means of full denoising, thus as closed-loop as the proposed method? Can the authors clarify this? 

The planning horizon for the two control tasks is fixed. Why exactly to these lengths? What happens if full denoising is done over shorter horizons, replanning at each time step, using DiffPhyCon-1?
soundness: 4
strengths: The main strength of the proposed method is the idea to save computation by not fully denoising the latter steps of the control sequence, with the understanding that in receding horizon control, the controls for those latter stages will likely not be used, as new state observations will appear in the meantime and the control sequence will be replanned. This idea addresses one of the main shortcomings of controllers based on diffusion policies - the very long process of denoising during sampling, which is currently a major obstacle to the wider adoption of diffusion policies that have shown considerable advantages with respect to more traditional behavior cloning schemes.
summary: The paper proposes a novel method for using diffusion policies for closed loop control that speeds up the generation of a control sequence of a length up to a chosen horizon by avoiding full denoising of the latter steps in the sequence, as they will likely not be executed by a receding horizon controller. The proposed method demonstrates superior performance in terms of achieving control objectives, due to closed loop control using a diffusion policy, with computation time much shorter than other diffusion policy baseline controllers.
weaknesses: I believe the proposed method is excellent and clearly applicable to very difficult control problems involving PDEs. However, most control problems are not of this kind. Why not show the performance on a simpler control problem with a handful of variables, hopefully not requiring high-end GPUs? This could potentially increase the applicability and appeal of the method significantly.

Some minor typos:

P. 3, L. 112: "It" -> "They"
P. 3, L. 140: "has" -> "have"
P. 7, L. 376: "well well" -> "perform well"?
P. 10, L. 539: "done" -> "drone"?
P. 19, L. 1009: "128 grids" -> "128 cells"

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 1
questions: A careful rewriting of the manuscript seems necessary to clarify the problem setting and contribution, and to appreciate the differences between the propose method and the existing ones. A characterization of the performance of the proposed method is necessary, otherwise the  method takes the form of a heuristic that may not be desirable in controls applications.

As a side note, the problem described in equation (1) is inherently an open loop problem, since no structure is imposed on the input w. In fact, the set of open-loop sequences include the set of closed-loop sequences, but the converse may not be true. Since the system dynamics are deterministic, it is also unclear why using a closed-loop strategy would be necessary. Finally, the used notation (w for control inputs and u for the system state) is very unusual and, to the best of my knowledge, never been used in any controls paper.
soundness: 2
strengths: The overall objective of the paper is interesting and timely: to design closed-loop diffusion-based control algorithms. The issue of computational complexity is important, since diffusion-based methods may work well but require extensive computations, making them difficult to use in feedback settings.
summary: The paper studies the problem of designing feedback control algorithms using diffusion techniques. Differently from open-loop control strategies, which compute a pre-determined control sequence that is applied to a system over time, closed-loop methods use information about the state of the system to tune the control inputs to compensate for deviations of the system trajectories from the predicted one, due to unmodelled dynamics or noise. The paper focuses on the computational complexity of the methods, and it proposes a technique to render closed-loop diffusion-based control computationally efficient.
weaknesses: The paper is poorly presented, to an extent that makes the contributions unclear and difficult to appreciate. This start with the definition of the problem (where what is known/unknown/available to the designer is unclear, where there is confusion between system and environment, where   the stated problem is inherently open-loop and the feedback structure is not evident), and continues into the description of the proposed method (the synchronous and asynchronous parts are vague and difficult to grasp). Additionally, and importantly, the paper doesn't seem to offer performance guarantees for the proposed method, which makes it difficult to advocate for the proposed method in any practical use.

Meta-review: The paper proposes a novel method for closed loop control using a diffusion policy which preferentially denoises early steps of the control sequence over later ones. The proposed method demonstrates superior performance in terms of achieving control objectives, with computation time much shorter than other diffusion policy baseline controllers.

The insight of the proposed method is simple (avoid denoising later steps of a controller since they are unlikely to be executed in a closed-loop system), but insightful. The result is an approach which is likely to be widely applicable to many control problems, as the authors demonstrate in several low and high-dimensional control settings. The proposed method achieves better performance results than prior work, while also having faster computation times.

The biggest weakness of the paper is in the presentation. The reviewers were split on the quality of the presentation, but in my own reading, the method *is* hard to follow. One reviewer also declined to provide a review at all because they could not understand the paper's contents despite reading it several times. Several variables are not defined, and Figure 2 (demonstrating the method) is not sufficiently self-encapsulating as an explanation on its own. For those already very familiar with diffusion models (and specifically how they might be used in a control setting), the math might be easier to follow, but it seems unlikely that a reader less well-versed in diffusion modeling (or control) could figure out how the method works without spending a week reading the paper in depth. In either case, the method is likely not described clearly enough for someone to reimplement it.

However, given that there were very positive reviews on the novelty and success of the method, I lean towards acceptance (poster). The authors are encouraged to take the concerns about presentation quality seriously for the camera-ready version however.

Decision: Accept (Poster)