Implementasi Praktis &
Integrasi SCADA / ICS

Bahasa Pemrograman IEC 61131-3 untuk Kontrol Adaptif

Kode berikut mensimulasikan konsep HIL: kontroller (dijalankan di "PLC virtual") dan plant simulator berjalan dalam loop yang sinkron, saling bertukar data melalui antarmuka yang meniru komunikasi OPC-UA. Ini adalah fondasi dari HIL testing sebelum deploy ke hardware nyata.

import numpy as np
import matplotlib.pyplot as plt
import time

# ==========================================
# SIMULASI HIL: GAIN SCHEDULING + MRAC
# untuk Kontrol Frekuensi Turbin Gas
# "Plant Simulator" ↔ "PLC Controller"
# ==========================================

# === PLANT SIMULATOR (mensimulasikan turbin gas) ===
class GasTurbinePlant:
    def __init__(self):
        self.freq  = 50.0     # frekuensi awal 50 Hz
        self.power = 0.0      # daya output MW (normalized)
        self.M     = 10.0     # inertia konstan
        self.D     = 1.5      # damping
        # Parameter berubah seiring waktu (simulasi degradasi)
        self.M_true = 10.0

    def step(self, u_gen, P_load, dt, t):
        # Parameter drift setelah t=100s (simulasi degradasi nyata)
        self.M_true = 10.0 - 0.02*(t > 100) * min(t-100, 150)
        df = (u_gen - P_load - self.D*(self.freq-50)) / self.M_true
        self.freq += df * dt
        return self.freq

# === PLC CONTROLLER (mensimulasikan kode Structured Text di PLC) ===
class AdaptivePLCController:
    def __init__(self):
        # Gain Scheduling: 3 zona berdasarkan |error frekuensi|
        self.zones = {
            'normal':  {'Kp':8.0,  'Ki':1.2, 'thresh':0.2},
            'moderate':{'Kp':14.0, 'Ki':2.0, 'thresh':0.6},
            'large':   {'Kp':22.0, 'Ki':3.5, 'thresh':9999},
        }
        self.integral = 0.0
        self.prev_e   = 0.0
        # MRAC layer untuk kompensasi drift
        self.theta    = 1.0    # adaptive gain multiplier
        self.gamma    = 0.005  # learning rate
        self.zone_hist = []

    def select_zone(self, err_abs):
        for name, z in self.zones.items():
            if err_abs <= z['thresh']: return name, z
        return 'large', self.zones['large']

    def compute(self, freq_meas, setpoint, dt):
        e      = setpoint - freq_meas
        e_abs  = abs(e)
        zone_name, zone = self.select_zone(e_abs)
        self.zone_hist.append(zone_name)

        # Gain Scheduling: Kp, Ki dari zona
        Kp_eff = zone['Kp'] * self.theta    # MRAC mengalikan gain
        Ki_eff = zone['Ki'] * self.theta

        self.integral = np.clip(self.integral + e*dt, -5, 5)
        u = np.clip(Kp_eff*e + Ki_eff*self.integral, -2.0, 2.0)

        # MRAC adaptation law: update theta berdasarkan tracking error
        # Jika |e| naik dari sebelumnya → theta perlu naik
        de = e - self.prev_e
        self.theta += self.gamma * e * de
        self.theta  = np.clip(self.theta, 0.5, 3.0)
        self.prev_e = e
        return u

# === SIMULASI LOOP HIL ===
dt       = 0.02
T_total  = 300
t_arr    = np.arange(0, T_total, dt)
N        = len(t_arr)

# Profil gangguan beban (MW, normalized)
P_load = np.zeros(N)
P_load[int(20/dt):int(60/dt)]   = 0.5   # load step +0.5
P_load[int(80/dt):int(120/dt)]  = -0.3  # load drop
P_load[int(150/dt):int(200/dt)] = 0.8   # besar setelah degradasi
P_load[int(230/dt):]              = 0.4

def run_hil(use_adaptive):
    plant  = GasTurbinePlant()
    ctrl   = AdaptivePLCController()
    f_hist, u_hist, theta_hist = [], [], []
    for k, t in enumerate(t_arr):
        f_meas = plant.freq + 0.005*np.random.randn()  # sensor noise
        if use_adaptive:
            u = ctrl.compute(f_meas, 50.0, dt)
        else:
            # PID statis baseline (satu zona saja)
            e = 50.0 - f_meas
            ctrl.integral = np.clip(ctrl.integral + e*dt, -5, 5)
            u = np.clip(10.0*e + 1.5*ctrl.integral, -2, 2)
        plant.step(u, P_load[k], dt, t)
        f_hist.append(plant.freq)
        u_hist.append(u)
        theta_hist.append(ctrl.theta)
    return np.array(f_hist), np.array(u_hist), np.array(theta_hist), ctrl.zone_hist

f_pid,   u_pid,   _, _            = run_hil(use_adaptive=False)
f_adp,   u_adp,   theta_h, zones  = run_hil(use_adaptive=True)

# === PLOT HASIL HIL ===
fig, axes = plt.subplots(3, 1, figsize=(13, 10))

axes[0].plot(t_arr, f_pid, color='tomato', lw=1.5, alpha=0.8, label='PID Statis')
axes[0].plot(t_arr, f_adp, color='lime',   lw=2.0,             label='GS + MRAC Adaptif')
axes[0].axhline(50, color='orange', ls='--', lw=1.5, label='Setpoint 50 Hz')
axes[0].fill_between(t_arr, 49.8, 50.2, alpha=0.07, color='lime', label='Zona ±0.2 Hz')
axes[0].axvline(100, color='yellow', ls=':', lw=1.5, alpha=0.6)
axes[0].annotate('Degradasi\ndimulai', xy=(100,50.3), color='yellow', fontsize=9)
axes[0].set_ylabel('Frekuensi (Hz)')
axes[0].set_title('HIL Test: GS+MRAC vs PID Statis — Kontrol Frekuensi Turbin Gas')
axes[0].legend(fontsize=9); axes[0].grid(alpha=0.3)

axes[1].plot(t_arr, theta_h, color='orange', lw=1.8, label='θ MRAC (adaptive gain multiplier)')
axes[1].axhline(1.0, color='white', ls='--', lw=1, alpha=0.4, label='θ=1 (nominal)')
axes[1].axvline(100, color='yellow', ls=':', lw=1.5, alpha=0.6)
axes[1].set_ylabel('θ (MRAC param)')
axes[1].set_title('θ Naik Otomatis setelah Degradasi: MRAC Mendeteksi & Mengkompensasi')
axes[1].legend(); axes[1].grid(alpha=0.3)

axes[2].plot(t_arr, P_load, color='red',    lw=2, label='Gangguan Beban')
axes[2].plot(t_arr, u_adp,  color='cyan',   lw=1.5, label='Output Kontrol Adaptif')
axes[2].plot(t_arr, u_pid,  color='tomato', lw=1.2, alpha=0.6, ls='--', label='Output PID Statis')
axes[2].set_xlabel('Waktu (s)'); axes[2].set_ylabel('Daya (p.u.)')
axes[2].set_title('Beban vs Output Kontrol — Adaptif lebih responsif di kondisi post-degradasi')
axes[2].legend(fontsize=9); axes[2].grid(alpha=0.3)

plt.tight_layout(); plt.show()

rmse_pid = np.sqrt(np.mean((f_pid - 50)**2))
rmse_adp = np.sqrt(np.mean((f_adp - 50)**2))
print(f"\n📊 RMSE Deviasi Frekuensi:")
print(f"   PID Statis : {rmse_pid:.4f} Hz")
print(f"   GS + MRAC  : {rmse_adp:.4f} Hz")
print(f"   Improvement: {(1-rmse_adp/rmse_pid)*100:.1f}%")
print(f"\n✅ θ final = {theta_h[-1]:.3f} — MRAC berhasil kompensasi degradasi turbin!")