4.1  Model-Free Methods

import numpy as np
import pandas as pd

np.random.seed(42)

# ══════════════════════════════════════════════════════════
# Q-Learning for Irrigation Scheduling
# State: soil_moisture (discretised 0–4), temp (0–2)
# Actions: 0=no water, 1=5mm, 2=10mm, 3=20mm
# Reward: yield gain - water cost
# ══════════════════════════════════════════════════════════
N_MOISTURE = 5   # 0–4 (very dry to saturated)
N_TEMP     = 3   # 0=cool, 1=warm, 2=hot
N_ACTIONS  = 4   # irrigation amounts
N_STATES   = N_MOISTURE * N_TEMP

# State encoding
def encode_state(moisture, temp):
    return moisture * N_TEMP + temp

# Simulate environment step
def step(moisture, temp, action):
    water_added = [0, 5, 10, 20][action]
    water_cost  = [0, 0.8, 1.5, 2.8][action]
    
    # Yield reward (peaks at moisture=3)
    new_moisture = min(4, max(0, moisture + action - int(0.3*temp+0.5)))
    yield_reward = max(0, 4 - abs(new_moisture - 3)) * 2.5
    reward       = yield_reward - water_cost
    
    # Next state transition
    new_temp = np.random.choice([0,1,2], p=[0.3,0.4,0.3])
    return encode_state(new_moisture, new_temp), reward, new_moisture, new_temp

# ── Initialise Q-Table ────────────────────────────────────
Q = np.zeros((N_STATES, N_ACTIONS))
alpha   = 0.1     # learning rate
gamma   = 0.95    # discount
epsilon = 1.0     # initial exploration
eps_min = 0.05
eps_decay = 0.999

episode_rewards = []

for episode in range(3000):
    moisture = np.random.randint(0, 5)
    temp     = np.random.randint(0, 3)
    state    = encode_state(moisture, temp)
    total_r  = 0
    
    for step_num in range(30):  # 30 days per episode
        # ε-greedy action selection
        if np.random.rand() < epsilon:
            action = np.random.randint(N_ACTIONS)
        else:
            action = np.argmax(Q[state])
        
        next_state, reward, moisture, temp = step(moisture, temp, action)
        
        # Q-Learning update (Bellman)
        td_error = reward + gamma * np.max(Q[next_state]) - Q[state, action]
        Q[state, action] += alpha * td_error
        
        state   = next_state
        total_r += reward
    
    epsilon = max(eps_min, epsilon * eps_decay)
    episode_rewards.append(total_r)

print("=== Q-Learning — Irrigation Scheduling ===")
print(f"Final ε = {epsilon:.4f}")
print(f"Last 200 episodes avg reward: {np.mean(episode_rewards[-200:]):.3f}")
print(f"First 200 episodes avg reward: {np.mean(episode_rewards[:200]):.3f}")

# ── Print optimal policy ─────────────────────────────────
print("\nOptimal Policy (state → action):")
action_names = ['No water','5mm','10mm','20mm']
print(f"{'Moisture':>10} | {'Temp':>6} | {'Best Action':>15} | {'Q-value':>10}")
for m in range(N_MOISTURE):
    for t in range(N_TEMP):
        s = encode_state(m, t)
        a = np.argmax(Q[s])
        temp_name = ['Cool','Warm','Hot'][t]
        print(f"{'Level '+str(m):>10} | {temp_name:>6} | {action_names[a]:>15} | {Q[s,a]:>10.4f}")

# ── Q-Table heat summary ─────────────────────────────────
print("\nQ-Table (max Q per state):")
print(np.round(Q.max(axis=1).reshape(N_MOISTURE, N_TEMP), 3))

import numpy as np

np.random.seed(42)

# ══════════════════════════════════════════════════════════
# REINFORCE for Greenhouse Control (simplified)
# State: [temp_error, co2_deficit] — 2 continuous dims
# Action: [delta_temp, delta_co2]  — 2 continuous
# Policy: Gaussian with learned mean, fixed std
# ══════════════════════════════════════════════════════════

class GaussianPolicy:
    def __init__(self, state_dim=2, action_dim=2, lr=0.01):
        np.random.seed(42)
        # Linear policy: mean = W @ state + b
        self.W = np.random.randn(action_dim, state_dim) * 0.1
        self.b = np.zeros(action_dim)
        self.std = 0.5
        self.lr  = lr

    def mean(self, state):
        return np.tanh(self.W @ state + self.b)  # bounded actions

    def sample(self, state):
        mu = self.mean(state)
        return np.random.normal(mu, self.std)

    def log_prob(self, state, action):
        """log pi_theta(a|s) up to a constant. Its gradient w.r.t. mu is
        (a - mu)/std**2 — exactly the `grad_mu` used in update() below,
        which is where the policy gradient theorem enters this code."""
        mu = self.mean(state)
        return -0.5 * np.sum(((action - mu) / self.std)**2)

    def update(self, states, actions, returns):
        dW = np.zeros_like(self.W)
        db = np.zeros_like(self.b)
        for s, a, G in zip(states, actions, returns):
            mu   = self.mean(s)
            # Policy gradient: ∇ log π × G
            grad_mu = (a - mu) / (self.std**2)
            # Chain through tanh
            dtanh = 1 - np.tanh(self.W @ s + self.b)**2
            dW += np.outer(grad_mu * dtanh, s) * G
            db += (grad_mu * dtanh) * G
        self.W += self.lr * dW / len(states)
        self.b += self.lr * db / len(states)

def greenhouse_env_step(state, action):
    """Simulated greenhouse: reward = growth - energy_cost"""
    temp_err  = state[0] - action[0] * 0.3 + np.random.normal(0, 0.1)
    co2_def   = state[1] - action[1] * 0.3 + np.random.normal(0, 0.1)
    growth    = max(0, 1.5 - temp_err**2 - 0.5*co2_def**2)
    energy    = 0.3 * np.sum(action**2)
    reward    = growth - energy
    next_state = np.clip([temp_err, co2_def], -2, 2)
    return next_state, reward

# ── Training REINFORCE ────────────────────────────────────
policy  = GaussianPolicy(lr=0.005)
gamma   = 0.99
ep_returns, ep_logp = [], []

for episode in range(1000):
    state = np.random.uniform(-1, 1, 2)
    states, actions, rewards = [], [], []
    
    for _ in range(20):  # 20 steps per episode
        action = policy.sample(state)
        next_state, reward = greenhouse_env_step(state, action)
        states.append(state.copy()); actions.append(action); rewards.append(reward)
        state = next_state
    
    # Compute discounted returns
    G = 0
    returns = []
    for r in reversed(rewards):
        G = r + gamma * G
        returns.insert(0, G)
    
    # Normalise returns (variance reduction)
    returns = np.array(returns)
    returns = (returns - returns.mean()) / (returns.std() + 1e-8)
    
    policy.update(states, actions, returns)
    ep_returns.append(sum(rewards))
    ep_logp.append(np.mean([policy.log_prob(st, ac)
                            for st, ac in zip(states, actions)]))

print("=== REINFORCE — Greenhouse Control ===")
print(f"Initial avg reward (ep 1-50):   {np.mean(ep_returns[:50]):.4f}")
print(f"Final   avg reward (ep 950-1000): {np.mean(ep_returns[950:]):.4f}")
print(f"Improvement: {100*(np.mean(ep_returns[950:])-np.mean(ep_returns[:50]))/abs(np.mean(ep_returns[:50])+1e-8):.1f}%")

# Mean log-probability of the actions actually taken: it rises as the policy
# concentrates mass on the actions it has learned are good.
print(f"Mean log pi(a|s)  — first 50: {np.mean(ep_logp[:50]):.4f}")
print(f"Mean log pi(a|s)  — last  50: {np.mean(ep_logp[-50:]):.4f}")

print(f"\nLearned policy weights W:\n{np.round(policy.W, 3)}")
print(f"Learned bias b: {np.round(policy.b, 3)}")