Surrogate-assisted DPC for a 2D Stefan phase-change material system under diurnal disturbances
Graduate research -- Johns Hopkins University, SOLARIS Lab. Advisors: Dr. Jan Drgona · Dr. Liang Wu. We train neural operator surrogates (DeepONet, FNO) to approximate the Stefan enthalpy PDE governing a PCM slab, then differentiate through them to solve the control problem via gradient descent -- following the DPC framework of Drgona et al. (2022).
| Controller | SoC RMSE | Pk-Pk SoC | Mean m_dot (g/s) | Re-plan (s) |
|---|---|---|---|---|
| Open-loop (constant flow) | 0.058 | 0.017 | 25.0 | < 0.001 |
| PI | 0.060 | 0.022 | 35.3 | < 0.001 |
| BC-Offline (neural policy) | 0.057 | 0.013 | 21.7 | 0.003 |
| DeepONet-DPC (proposed) | 0.053 | 0.004 | 4.9 | 1.5 |
| FNO-DPC | 0.053 | 0.009 | 28.0 | 5.0 |
| MPPI-FNO (K=200) | 0.051 | 0.004 | 17.6 | 0.8 |
| PhiFlow-DPC (H=5) | 0.053 | 0.005 | 8.8 | 5.0 |
| PhiFlow-DPC (H=10, ceiling) | 0.050 | 0.003 | 6.0 | 13.1 |
DeepONet-DPC matches PhiFlow-DPC at matched horizon (0.053 vs 0.053) at 3.3x faster re-plan time, and is the only gradient-based controller that fits within a 30 s real-time control interval.
PCM thermal energy storage systems absorb and release heat through a solid-liquid phase transition, governed by the Stefan enthalpy PDE. The phase change creates a sharp nonlinearity (Gaussian peak in effective heat capacity) that breaks linear MPC. Full PDE solvers are too slow for real-time optimization (30-300 s per step in 2D).
This project:
- Builds a FEM-validated 2D finite-difference simulator for Rubitherm RT37HC paraffin (15x15 grid, explicit Euler, dt=30 s)
- Trains DeepONet and FNO surrogates to approximate the PDE in a single forward pass (~1 ms)
- Implements DPC by backpropagating a receding-horizon SoC tracking cost through the surrogate
- Benchmarks 9 controllers with harmonized cost functions on a 2-hour regulation scenario
- Trains a behaviour-cloned offline policy (~2 ms/decision) from online DPC demonstrations
| Parameter | Value |
|---|---|
| PCM material | Rubitherm RT37HC paraffin |
| Melting temperature | 310.15 K (37 C) |
| Latent heat | 210 kJ/kg |
| Domain | 0.1 x 0.1 m, 15x15 grid |
| Control input | HTF mass flow rate, 0.001-0.05 kg/s |
| Disturbance | Sinusoidal inlet temperature + AR(1) noise |
| Control objective | Regulate SoC = 0.50 |
| Time step | 30 s |
The HTF flows along the left wall with convective BC; the other three walls are insulated. The control problem is a receding-horizon optimization:
min sum_{h=0}^{H-1} Q*(SoC_h - 0.5)^2 + R*m_dot_h^2 + R2*(m_dot_h - m_dot_{h-1})^2
The NumPy FDM simulator was validated through a four-level cascade before use as a data generator:
| Level | Test | Result |
|---|---|---|
| 1 | Analytical BC energy balance | < 1e-10 K error |
| 2 | Cross-simulator (NumPy vs PhiFlow) | 0.001 K (pure liquid), 0.133 K (mushy) |
| 3 | SoC enthalpy self-consistency | < 0.001 |
| 4 | FEM ground-truth (FEniCS 40x40, Crank-Nicolson) | SoC RMSE = 0.0008 |
git clone https://github.com/SOLARIS-JHU/PCM-DPC.git
cd PCM-DPC
pip install -r requirements.txtGPU (CUDA) recommended. Tested on RTX 4070 Laptop (CUDA 12).
python run_validation.pypython -c "
from pcm_datagen.dataset_2d import generate_dataset_2d
generate_dataset_2d(N=1000, output_dir='output/dataset_2d', n_workers=6, consolidate=True)
"python src/train_deeponet_2d.py --device cuda --epochs 200 --lr 3e-4 \
--data_dir output/dataset_2d --save models/deeponet_2d.ptpython src/run_comparison_5way_2d.py --device cuda# Collect DPC decisions (~46 min, 150 episodes)
python src/collect_bc_data.py --n_episodes 150 --H 5 --n_iter 20 --seed 42
# Train via behaviour cloning (~25 sec, 200 epochs)
python src/train_bc_policy.py --epochs 200python src/run_r2_sweep.py --r2_values 0 0.01 0.05 0.1 0.2 0.5 1.0 2.0
python src/run_horizon_sweep.py --device cuda --H_values 5 10 20 40
python src/run_diurnal_comparison.py --device cudaPCM-DPC/
├── src/ # Controllers and training
│ ├── dpc_controller_2d.py # Online DPC via DeepONet
│ ├── dpc_controller_2d_fno.py # Online DPC via FNO
│ ├── dpc_controller_2d_phiflow.py # DPC via differentiable PhiFlow PDE
│ ├── dpc_controller_2d_mppi.py # MPPI sampling-based controller
│ ├── dpc_controller_2d_offline.py # BC-Offline policy inference
│ ├── train_deeponet_2d.py # DeepONet surrogate training
│ ├── train_bc_policy.py # Behaviour cloning training
│ ├── collect_bc_data.py # Collect (state -> m_opt) pairs from online DPC
│ ├── neural_operator_2d.py # DeepONet2D architecture
│ ├── neural_operator_benchmarks.py # FNO/CNO architectures
│ ├── run_comparison_5way_2d.py # 9-controller benchmark
│ ├── run_r2_sweep.py # R2 smoothness penalty sweep
│ ├── run_horizon_sweep.py # Prediction horizon H sweep
│ └── run_diurnal_comparison.py # 24-hour diurnal comparison
│
├── pcm_datagen/ # Data generation (FEM-validated)
│ ├── simulator_2d.py # Stefan PDE solver (FDM, NumPy)
│ ├── simulator_2d_phiflow.py # Differentiable PhiFlow solver
│ ├── dataset_2d.py # Trajectory generator
│ └── signals.py # Disturbance signal generators
│
├── fem_benchmark/
│ └── fenics_pcm_2d.py # FEniCS FEM oracle (validation only)
│
├── scripts/ # Post-processing
│ ├── analyze_comparison.py
│ ├── make_animation_2d.py
│ └── make_animation_24h.py
│
├── config.py # Canonical experiment parameters
├── run_validation.py # 4-level validation cascade
├── requirements.txt
└── docs/ # Paper and validation reports
models/andoutput/are git-ignored (large binaries). Pre-trained weights and datasets available on request.
The surrogate maps the current temperature field + control input to the next temperature field in one forward pass:
T_hat_{t+1}(x) = branch(T_t, m_dot, T_in) . trunk(x)^T + b
- Branch net: MLP [227, 256, 256, 128, 128], ReLU
- Trunk net: MLP [3, 64, 128, 128, 128], Tanh
- Trunk augmentation: input is [x, y, exp(-x/delta)] where delta = sqrt(alpha*dt) is the thermal diffusion length -- this encodes the boundary-layer structure at the convective wall
- Parameters: ~440,000
- Training: 1,000 trajectories, 1.38M state-transition pairs, Adam lr=3e-4, 200 epochs
The control sequence u = [m_dot_0, ..., m_dot_{H-1}] is treated as a free parameter. Automatic differentiation computes dJ/du through the unrolled surrogate, and Adam updates u at each re-plan step. SoC is computed with a soft-sigmoid (hard clip zeros the gradient at saturation). Input constraints are enforced via sigmoid reparameterization.
- PhiFlow-DPC: Differentiates through real PDE steps (no surrogate error). Performance ceiling.
- DeepONet-DPC: Differentiates through the surrogate. 3.3x faster than PhiFlow-DPC.
- MPPI-FNO: Sampling-based (K=200 rollouts through FNO). No gradients needed.
- FNO-DPC: Gradient-based DPC through FNO surrogate.
- BC-Offline: Single forward pass (~2 ms). Trained from DeepONet-DPC demonstrations.
- PI, Bang-Bang, Open-loop: Classical baselines.
Level 4 validation requires FEniCS via Docker:
docker pull quay.io/fenicsproject/stable
docker run -v $(pwd):/home/fenics/shared quay.io/fenicsproject/stable \
python /home/fenics/shared/fem_benchmark/fenics_pcm_2d.py@article{drgona2022dpc,
title={Differentiable Predictive Control: An MPC Alternative for Unknown Nonlinear Systems using Constrained Deep Learning},
author={Drgo{\v{n}}a, J{\'a}n and Kis, Karol and Tuor, Aaron and Vrabie, Draguna and Klauco, Martin},
journal={Journal of Process Control},
volume={116},
pages={80--92},
year={2022}
}
@article{lu2021deeponet,
title={Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators},
author={Lu, Lu and Jin, Pengzhan and Pang, Guofei and Zhang, Zhongqiang and Karniadakis, George Em},
journal={Nature Machine Intelligence},
volume={3},
pages={218--229},
year={2021}
}
@article{voller1990stefan,
title={Fixed grid techniques for phase change problems: a review},
author={Voller, Vaughan R},
journal={International Journal for Numerical Methods in Engineering},
year={1990}
}Tools: PyTorch · PhiFlow · neuraloperator · FEniCS · NumPy
Jose Maria Borrego Acosta -- Graduate student, Johns Hopkins University Jan Drgona -- Johns Hopkins University Liang Wu -- Johns Hopkins University
MIT