# PI(D)-regeling van de stationtemperatuur: simulatie met de methode van Euler
TAU = 10800.0     # tijdconstante van het station (s): 3 uur
K = 0.02          # statische versterking (graden Celsius per watt)
U_MAX = 2000.0    # vermogen van de verwarming (W)


def simuleer(Kp, Ki, Kd=0.0, r=15.0, T=32400, dt=10.0, anti_windup=True):
    y = 0.0               # temperatuurstijging boven de buitentemperatuur (graden Celsius)
    integraal = 0.0
    vorige = y
    tijden, temps = [], []
    for k in range(int(T / dt) + 1):
        e = r - y
        afgeleide = -(y - vorige) / dt           # D op de meting: geen sprong bij een nieuw setpoint
        u_vrij = Kp * e + Ki * integraal + Kd * afgeleide
        u = min(max(u_vrij, 0.0), U_MAX)          # de verwarming kan enkel 0 tot 2000 W leveren
        verzadigd = u != u_vrij
        if not (anti_windup and verzadigd and (u_vrij > 0) == (e > 0)):
            integraal += e * dt                   # anti-windup: niet verder integreren tegen de grens
        vorige = y
        y += dt * (-y + K * u) / TAU              # Euler-stap van tau * dy/dt = -y + K * u
        tijden.append(k * dt)
        temps.append(y)
    return tijden, temps


t, y = simuleer(Kp=200, Ki=0.17)
print("met anti-windup:   maximum %.2f, na negen uur %.2f" % (max(y), y[-1]))
t, y = simuleer(Kp=200, Ki=0.17, anti_windup=False)
print("zonder anti-windup: maximum %.2f, na negen uur %.2f" % (max(y), y[-1]))