""" Reinforcement-Learning-Augmented Control Simulation for a PV-Integrated Welding Power Supply ---------------------------------------------------- This is a NEW, simplified simulation built to define a concrete, reproducible RL methodology (algorithm, state space, action space, reward function) for the paper. It is NOT a reconstruction of the original hardware experiment reported in the source chapter. Implemented from scratch with numpy only (no torch/gym available, offline environment). Architecture matches the paper's Section 4 description: - Inner loop: fast PI controller regulating arc current (represented here as a first-order closed-loop tracking model with PI gains). - Outer loop: RL agent (DQN) that outputs a small discrete trim to the inner loop's current reference every outer step, in order to minimize a composite cost (tracking error, DC-link ripple, efficiency loss proxy, THD proxy) under PV intermittency and pulsed welding load disturbances. """ import numpy as np import json rng = np.random.default_rng(42) # ---------------------------------------------------------------------- # 1. ENVIRONMENT # ---------------------------------------------------------------------- class WeldingPVEnv: """ Simplified discrete-time model of the arc-current / DC-link dynamics of a PV-fed welding power converter, driven by: - an irradiance profile (PV available power fraction, 0-1) - a pulsed welding load profile (reference arc current, p.u.) State (normalized, in [-1, 1] roughly): [ i_err, v_ripple, pv_ratio, load_state ] Action (discrete, 5 choices): outer-loop trim added to the inner PI current reference {-0.02, -0.01, 0.0, +0.01, +0.02} p.u. """ ACTIONS = np.array([0.4, 0.7, 1.0, 1.5, 2.2]) # gain-scheduling multipliers N_ACTIONS = len(ACTIONS) STATE_DIM = 4 KP0, KI0 = 0.9, 0.04 NOMINAL_ARC_LEN = 0.65 NOMINAL_GAIN = 0.12 + 0.45 * (1 - NOMINAL_ARC_LEN) # process gain at nominal arc length def __init__(self, episode_len=400, dt=1.0, seed=None): self.episode_len = episode_len self.dt = dt self.rng = np.random.default_rng(seed) self.reset() def _irradiance_profile(self, t): # slow-varying cloud-induced fluctuation + fast noise base = 0.75 + 0.20 * np.sin(2 * np.pi * t / 250.0) cloud_dip = 0.25 * (np.sin(2 * np.pi * t / 90.0) < -0.85) noise = self.rng.normal(0, 0.02) return float(np.clip(base - cloud_dip + noise, 0.05, 1.0)) def _load_profile(self, t): # pulsed welding current reference (p.u. of rated arc current) pulse = 0.5 + 0.5 * (np.sin(2 * np.pi * t / 40.0) > 0.0) ramp_noise = self.rng.normal(0, 0.01) return float(np.clip(pulse + ramp_noise, 0.1, 1.0)) def _update_arc_length(self): # slow mean-reverting drift representing torch-handling variation, # plus occasional abrupt operator-induced jumps self.arc_len += 0.06 * (self.NOMINAL_ARC_LEN - self.arc_len) + self.rng.normal(0, 0.02) if self.rng.random() < 0.012: self.arc_len = self.rng.uniform(0.25, 1.0) self.arc_len = float(np.clip(self.arc_len, 0.25, 1.0)) def reset(self): self.t = 0 self.i_arc = 0.5 # actual arc current, p.u. self.v_link = 1.0 # DC-link voltage, p.u. self.integral_err = 0.0 self.arc_len = self.NOMINAL_ARC_LEN self.pv_ratio = self._irradiance_profile(self.t) self.i_ref = self._load_profile(self.t) return self._get_state() def _get_state(self): i_err = self.i_ref - self.i_arc v_ripple = self.v_link - 1.0 arc_len_meas = float(np.clip(self.arc_len + self.rng.normal(0, 0.05), 0.0, 1.2)) return np.array([i_err, v_ripple, arc_len_meas - 0.65, self.pv_ratio - 0.5], dtype=np.float64) def step(self, action_idx, supply_mode="direct"): gain_mult = self.ACTIONS[action_idx] self._update_arc_length() process_gain = 0.12 + 0.45 * (1.0 - self.arc_len) # time-varying process gain # --- outer-loop-scheduled inner PI loop --- Kp, Ki = self.KP0 * gain_mult, self.KI0 * gain_mult err = self.i_ref - self.i_arc self.integral_err = np.clip(self.integral_err + err * self.dt, -0.6, 0.6) control = Kp * err + Ki * self.integral_err # supply behaviour: direct PV-only vs battery-buffered PV if supply_mode == "battery": if not hasattr(self, "_batt_state"): self._batt_state = self.pv_ratio self._batt_state += 0.15 * (self.pv_ratio - self._batt_state) effective_pv = self._batt_state batt_loss_penalty = 0.05 else: effective_pv = self.pv_ratio batt_loss_penalty = 0.0 pv_support = np.clip(effective_pv / 0.75, 0.3, 1.3) self.i_arc += self.dt * process_gain * pv_support * control self.i_arc = float(np.clip(self.i_arc, 0.0, 1.3)) # DC-link voltage ripple grows with mismatch between generation # and demand, and with abrupt trims (extra switching stress) mismatch = (self.i_arc - effective_pv) self.v_link += self.dt * (0.05 * mismatch - 0.10 * (self.v_link - 1.0)) self.v_link += self.rng.normal(0, 0.003) # advance disturbances self.t += self.dt self.pv_ratio = self._irradiance_profile(self.t) self.i_ref = self._load_profile(self.t) # --- reward: weighted composite cost (to be minimized) --- track_err = abs(self.i_ref - self.i_arc) ripple = abs(self.v_link - 1.0) gain_dev = abs(gain_mult - 1.0) thd_proxy = 0.5 * abs(control) + 0.3 * gain_dev # switching harshness eff_loss_proxy = 0.10 * gain_dev + 0.3 * ripple + 0.2 * track_err + batt_loss_penalty w = dict(track=4.0, ripple=2.0, thd=1.0, eff=1.5) cost = (w["track"] * track_err + w["ripple"] * ripple + w["thd"] * thd_proxy + w["eff"] * eff_loss_proxy) reward = -cost * 0.5 done = self.t >= self.episode_len info = dict(track_err=track_err, ripple=ripple, thd_proxy=thd_proxy, eff_loss_proxy=eff_loss_proxy) return self._get_state(), reward, done, info # ---------------------------------------------------------------------- # 2. DQN AGENT (small MLP, manual forward/backward with numpy) # ---------------------------------------------------------------------- class MLP: def __init__(self, sizes, seed=0): r = np.random.default_rng(seed) self.W, self.b = [], [] for i in range(len(sizes) - 1): limit = np.sqrt(6 / (sizes[i] + sizes[i + 1])) self.W.append(r.uniform(-limit, limit, (sizes[i], sizes[i + 1]))) self.b.append(np.zeros(sizes[i + 1])) def forward(self, x): acts = [x] for i in range(len(self.W) - 1): x = np.maximum(0, x @ self.W[i] + self.b[i]) # ReLU acts.append(x) x = x @ self.W[-1] + self.b[-1] # linear output acts.append(x) return x, acts def backward(self, acts, dY, lr, clip=5.0): dY = np.clip(dY, -clip, clip) grad = dY for i in reversed(range(len(self.W))): a_prev = acts[i] dW = np.clip(a_prev.T @ grad / grad.shape[0], -clip, clip) db = np.clip(grad.mean(axis=0), -clip, clip) self.W[i] -= lr * dW self.b[i] -= lr * db if i > 0: grad = grad @ self.W[i].T grad = grad * (acts[i] > 0) grad = np.clip(grad, -clip, clip) def copy_from(self, other, tau=1.0): for i in range(len(self.W)): self.W[i] = tau * other.W[i] + (1 - tau) * self.W[i] self.b[i] = tau * other.b[i] + (1 - tau) * self.b[i] class DQNAgent: def __init__(self, state_dim, n_actions, seed=0): hidden = 24 self.q = MLP([state_dim, hidden, hidden, n_actions], seed=seed) self.q_target = MLP([state_dim, hidden, hidden, n_actions], seed=seed + 1) self.q_target.copy_from(self.q, tau=1.0) self.n_actions = n_actions self.buffer = [] self.buffer_cap = 20000 self.gamma = 0.95 self.lr = 0.003 self.batch = 64 def act(self, state, eps): if rng.random() < eps: return rng.integers(self.n_actions) q_vals, _ = self.q.forward(state[None, :]) return int(np.argmax(q_vals[0])) def remember(self, s, a, r, s2, done): if len(self.buffer) >= self.buffer_cap: self.buffer.pop(0) self.buffer.append((s, a, r, s2, done)) def train_step(self): if len(self.buffer) < self.batch: return idx = rng.choice(len(self.buffer), self.batch, replace=False) S = np.array([self.buffer[i][0] for i in idx]) A = np.array([self.buffer[i][1] for i in idx]) R = np.array([self.buffer[i][2] for i in idx]) S2 = np.array([self.buffer[i][3] for i in idx]) D = np.array([self.buffer[i][4] for i in idx], dtype=np.float64) q_next, _ = self.q_target.forward(S2) target = R + self.gamma * (1 - D) * np.max(q_next, axis=1) q_vals, acts = self.q.forward(S) dY = np.zeros_like(q_vals) pred = q_vals[np.arange(self.batch), A] error = (pred - target) dY[np.arange(self.batch), A] = error self.q.backward(acts, dY, self.lr) def update_target(self, tau=0.05): self.q_target.copy_from(self.q, tau=tau) # ---------------------------------------------------------------------- # 3. TRAINING # ---------------------------------------------------------------------- def train(n_episodes=400): env = WeldingPVEnv(episode_len=400, seed=1) agent = DQNAgent(env.STATE_DIM, env.N_ACTIONS, seed=0) rewards_hist = [] eps = 1.0 for ep in range(n_episodes): s = env.reset() ep_reward = 0.0 done = False while not done: a = agent.act(s, eps) s2, r, done, info = env.step(a) agent.remember(s, a, r, s2, done) agent.train_step() s = s2 ep_reward += r agent.update_target(tau=0.05) eps = max(0.05, eps * 0.99) rewards_hist.append(ep_reward) if (ep + 1) % 25 == 0: print(f"Episode {ep+1:4d} | reward {ep_reward:8.2f} | eps {eps:.3f}") return env, agent, rewards_hist # ---------------------------------------------------------------------- # 4. FOUR CONTROL STRATEGIES # ---------------------------------------------------------------------- def rule_based_action(state): """Simple deterministic gain-scheduling heuristic (no learning), based only on the measured arc length (state[2], offset from nominal 0.65). Short arc -> high process gain -> schedule LOWER controller gain. Long arc -> low process gain -> schedule HIGHER controller gain.""" arc_len_dev = state[2] # arc_len_meas - 0.65 if arc_len_dev < -0.30: return 0 # mult = 0.4 (short arc, reduce gain) elif arc_len_dev < -0.10: return 1 # mult = 0.7 elif arc_len_dev < 0.10: return 2 # mult = 1.0 elif arc_len_dev < 0.25: return 3 # mult = 1.5 return 4 # mult = 2.2 (long arc, boost gain) def run_episode(env_seed, method, agent=None): """method in {'PI', 'Rule', 'Battery', 'RL'}""" supply_mode = "battery" if method == "Battery" else "direct" env = WeldingPVEnv(episode_len=400, seed=env_seed) s = env.reset() log = dict(t=[], i_ref=[], i_arc=[], v_link=[], track_err=[], ripple=[], thd_proxy=[], eff_loss_proxy=[]) done = False while not done: if method == "PI" or method == "Battery": a = 2 # trim = 0, inner PI only elif method == "Rule": a = rule_based_action(s) elif method == "RL": a = agent.act(s, eps=0.0) else: raise ValueError(method) log["t"].append(env.t) log["i_ref"].append(env.i_ref) s, r, done, info = env.step(a, supply_mode=supply_mode) log["i_arc"].append(env.i_arc) log["v_link"].append(env.v_link) for k in ("track_err", "ripple", "thd_proxy", "eff_loss_proxy"): log[k].append(info[k]) return log def evaluate(env_seed, method, agent=None): log = run_episode(env_seed, method, agent) return {k: float(np.mean(log[k])) for k in ("track_err", "ripple", "thd_proxy", "eff_loss_proxy")}, log def eff_pct(eff_loss_proxy): return 100 * (1 - eff_loss_proxy) if __name__ == "__main__": env, agent, rewards_hist = train(n_episodes=400) methods = ["Battery", "PI", "Rule", "RL"] results, logs = {}, {} for m in methods: metrics, log = evaluate(env_seed=999, method=m, agent=agent) logs[m] = log results[m] = { "mean_abs_tracking_error_pu": metrics["track_err"], "mean_abs_dc_link_ripple_pu": metrics["ripple"], "thd_proxy": metrics["thd_proxy"], "efficiency_proxy_pct": eff_pct(metrics["eff_loss_proxy"]), } print("\n=== Evaluation over fixed disturbance realization (seed=999) ===") print(json.dumps(results, indent=2)) with open("/home/claude/rl_results.json", "w") as f: json.dump({"training_rewards": rewards_hist, "evaluation": results}, f, indent=2) # ------------------------------------------------------------------ # FIGURES # ------------------------------------------------------------------ import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt colors = {"Battery": "#8c8c8c", "PI": "#4c72b0", "Rule": "#dd8452", "RL": "#55a868"} # --- Figure 1: bar comparison across 4 metrics --- fig, axes = plt.subplots(1, 4, figsize=(16, 4)) metric_keys = ["mean_abs_tracking_error_pu", "mean_abs_dc_link_ripple_pu", "thd_proxy", "efficiency_proxy_pct"] titles = ["Mean |Tracking Error| (p.u.)", "Mean |DC-link Ripple| (p.u.)", "THD proxy (a.u.)", "Efficiency proxy (%)"] for ax, key, title in zip(axes, metric_keys, titles): vals = [results[m][key] for m in methods] bars = ax.bar(methods, vals, color=[colors[m] for m in methods]) ax.set_title(title, fontsize=10) ax.tick_params(axis='x', rotation=20) for b, v in zip(bars, vals): ax.text(b.get_x() + b.get_width() / 2, v, f"{v:.3f}", ha="center", va="bottom", fontsize=8) fig.suptitle("Simulated Comparison: Battery vs PI vs Rule-Based vs RL-Augmented Control", fontsize=12) fig.tight_layout() fig.savefig("/home/claude/fig_comparison_bars.png", dpi=200) plt.close(fig) # --- Figure 2: time-series tracking comparison (first 150 steps) --- fig, axes = plt.subplots(4, 1, figsize=(10, 10), sharex=True) window = 150 for ax, m in zip(axes, methods): log = logs[m] t = log["t"][:window] ax.plot(t, log["i_ref"][:window], "k--", lw=1, label="Reference $i_{ref}$") ax.plot(t, log["i_arc"][:window], color=colors[m], lw=1.5, label=f"Arc current ({m})") ax.set_ylabel("Current (p.u.)") ax.legend(loc="upper right", fontsize=8) ax.set_title(m, fontsize=10, loc="left") axes[-1].set_xlabel("Time step") fig.suptitle("Arc Current Tracking Under Identical PV/Load Disturbance", fontsize=12) fig.tight_layout() fig.savefig("/home/claude/fig_tracking_timeseries.png", dpi=200) plt.close(fig) # --- Figure 3: training reward curve --- fig, ax = plt.subplots(figsize=(7, 4)) ax.plot(rewards_hist, color="#55a868") ax.set_xlabel("Training episode") ax.set_ylabel("Episode cumulative reward") ax.set_title("DQN Training Convergence (Outer-Loop RL Policy)") fig.tight_layout() fig.savefig("/home/claude/fig_training_curve.png", dpi=200) plt.close(fig) print("\nSaved figures: fig_comparison_bars.png, fig_tracking_timeseries.png, fig_training_curve.png")