Neural Network untuk
Identifikasi & Kontrol Adaptif

Teorema Universal Approximation (Cybenko, 1989) menyatakan bahwa satu hidden layer NN dengan neuron yang cukup dapat mengaproksimasi fungsi kontinu apapun dengan akurasi sembarang. Ini dasar mengapa NN sangat kuat untuk identifikasi sistem nonlinier.

import numpy as np
import matplotlib.pyplot as plt

# ==========================================
# NN IDENTIFIER DARI NOL (hanya NumPy)
# Identifikasi dinamika turbin gas nonlinier
# ==========================================

np.random.seed(42)

# === GENERATE DATA SISTEM NYATA (nonlinier) ===
# Turbin gas: y(k+1) = 0.7*y(k) + (1-exp(-u(k))) + 0.1*y(k)*u(k) + noise
def true_system(y_prev, u):
    return 0.7*y_prev + (1 - np.exp(-u)) + 0.1*y_prev*u

N_train = 800; N_test = 200; N = N_train + N_test
u_all = 2*np.random.rand(N) - 0.5   # input bervariasi (-0.5 to 1.5)
y_all = np.zeros(N); y_all[0] = 0.5
for k in range(1, N):
    y_all[k] = true_system(y_all[k-1], u_all[k-1]) + 0.02*np.random.randn()

# === BUAT DATASET: Input = [y(k), u(k)], Output = y(k+1) ===
X = np.column_stack([y_all[:-1], u_all[:-1]])
Y = y_all[1:]
X_train, Y_train = X[:N_train], Y[:N_train]
X_test,  Y_test  = X[N_train:], Y[N_train:]

# === NORMALISASI ===
X_mean, X_std = X_train.mean(0), X_train.std(0) + 1e-8
Y_mean, Y_std = Y_train.mean(), Y_train.std() + 1e-8
Xn_tr = (X_train - X_mean) / X_std
Yn_tr = (Y_train - Y_mean) / Y_std
Xn_te = (X_test  - X_mean) / X_std

# === NEURAL NETWORK: 2 → 16 → 1 ===
n_in=2; n_hid=16; n_out=1
W1 = np.random.randn(n_in, n_hid) * 0.3
b1 = np.zeros(n_hid)
W2 = np.random.randn(n_hid, n_out) * 0.3
b2 = np.zeros(n_out)

def relu(z): return np.maximum(0, z)
def relu_d(z): return (z > 0).astype(float)

def forward(X):
    z1 = X @ W1 + b1; h1 = relu(z1)
    z2 = h1 @ W2 + b2
    return z2.flatten(), z1, h1

# === TRAINING (Batch Gradient Descent) ===
lr = 0.005; epochs = 300; losses = []
for ep in range(epochs):
    yp, z1, h1 = forward(Xn_tr)
    loss = np.mean((yp - Yn_tr)**2)
    losses.append(loss)
    
    # Backprop
    dL_dyp = 2*(yp - Yn_tr) / len(Yn_tr)
    dL_dW2 = h1.T @ dL_dyp.reshape(-1,1)
    dL_db2 = dL_dyp.sum()
    dL_dh1 = dL_dyp.reshape(-1,1) @ W2.T
    dL_dz1 = dL_dh1 * relu_d(z1)
    dL_dW1 = Xn_tr.T @ dL_dz1
    dL_db1 = dL_dz1.sum(0)
    
    W1 -= lr*dL_dW1; b1 -= lr*dL_db1
    W2 -= lr*dL_dW2; b2 -= lr*dL_db2
    
    if ep % 50 == 0:
        print(f"Epoch {ep:3d}: Loss = {loss:.6f}")

# === TEST ===
yp_test_n, _, _ = forward(Xn_te)
yp_test = yp_test_n * Y_std + Y_mean  # denormalisasi
rmse = np.sqrt(np.mean((yp_test - Y_test)**2))

# === PLOT ===
fig, axes = plt.subplots(3, 1, figsize=(12, 9))

axes[0].semilogy(losses, color='orange', lw=1.8)
axes[0].set_xlabel('Epoch'); axes[0].set_ylabel('MSE Loss (log scale)')
axes[0].set_title('Training Loss — NN Belajar Dinamika Turbin Gas')
axes[0].grid(alpha=0.3)

t_test = np.arange(len(Y_test))
axes[1].plot(t_test, Y_test,  'b',      lw=2,   label='Output Aktual (sistem nyata)')
axes[1].plot(t_test, yp_test, 'orange', lw=1.5, ls='--', label=f'Prediksi NN (RMSE={rmse:.4f})')
axes[1].set_ylabel('Output y')
axes[1].set_title('Test Set: Output Aktual vs Prediksi NN')
axes[1].legend(); axes[1].grid(alpha=0.3)

err_test = Y_test - yp_test
axes[2].plot(t_test, err_test, 'tomato', lw=1.2, label='Error prediksi')
axes[2].axhline(0, color='lime', ls='--')
axes[2].fill_between(t_test, err_test, 0, alpha=0.2, color='tomato')
axes[2].set_xlabel('Waktu (k)'); axes[2].set_ylabel('Error')
axes[2].set_title('Residual Error — Mendekati Nol = NN Akurat!')
axes[2].legend(); axes[2].grid(alpha=0.3)

plt.tight_layout(); plt.show()
print(f"\n✅ RMSE Test: {rmse:.4f}")
print(f"NN berhasil mempelajari dinamika nonlinier turbin tanpa model fisika!")