# Digital Control

Biquad-cascade IIR filters, PID controller, lead/lag compensator design, and PID tuning rules. All components are no-std compatible, use no complex arithmetic at runtime, and work with both `f32` and `f64`.

Requires the `control` Cargo feature:

```toml
numeris = { version = "0.5", features = ["control"] }
```

## IIR Filter Design

numeris designs IIR filters by:

1. Computing analog prototype poles (Butterworth or Chebyshev Type I)
2. Applying a bilinear transform (Tustin's method) with frequency pre-warping
3. Pairing poles into second-order sections (biquads)

### Butterworth

Maximally flat passband, monotone rolloff.

```rust
use numeris::control::{butterworth_lowpass, butterworth_highpass, BiquadCascade};

// 4th-order lowpass at 1 kHz, 8 kHz sample rate
// → 2 biquad sections (N/2 = 2 for even order)
let mut lpf: BiquadCascade<f64, 2> = butterworth_lowpass(4, 1000.0, 8000.0).unwrap();

// Process one sample
let y = lpf.tick(1.0);

// Process a block
let input  = [1.0_f64, 0.5, -0.3, 0.8, 0.2];
let mut output = [0.0_f64; 5];
lpf.reset();
lpf.process(&input, &mut output);

// In-place processing
let mut buf = [1.0_f64, 0.5, -0.3];
lpf.reset();
lpf.process_inplace(&mut buf);

// Highpass version
let mut hpf: BiquadCascade<f64, 2> = butterworth_highpass(4, 2000.0, 8000.0).unwrap();
```

### Chebyshev Type I

Steeper rolloff than Butterworth at the cost of passband ripple.

```rust
use numeris::control::{chebyshev1_lowpass, chebyshev1_highpass, BiquadCascade};

// 4th-order Chebyshev lowpass, 1 dB passband ripple
let mut cheb: BiquadCascade<f64, 2> = chebyshev1_lowpass(4, 1.0, 1000.0, 8000.0).unwrap();

let y = cheb.tick(1.0);

// Highpass version
let mut cheb_hp: BiquadCascade<f64, 2> = chebyshev1_highpass(4, 1.0, 2000.0, 8000.0).unwrap();
```

### Design Function Signatures

```rust
// Butterworth
fn butterworth_lowpass<T, const N: usize>(
    order: usize,           // filter order (1..=2N, even fills N sections, odd uses degenerate last)
    cutoff_hz: T,           // -3 dB frequency
    sample_hz: T,           // sample rate
) -> Result<BiquadCascade<T, N>, ControlError>

fn butterworth_highpass<T, const N: usize>(
    order: usize,
    cutoff_hz: T,
    sample_hz: T,
) -> Result<BiquadCascade<T, N>, ControlError>

// Chebyshev Type I
fn chebyshev1_lowpass<T, const N: usize>(
    order: usize,
    ripple_db: T,           // passband ripple in dB (> 0)
    cutoff_hz: T,
    sample_hz: T,
) -> Result<BiquadCascade<T, N>, ControlError>

fn chebyshev1_highpass<T, const N: usize>(
    order: usize,
    ripple_db: T,
    cutoff_hz: T,
    sample_hz: T,
) -> Result<BiquadCascade<T, N>, ControlError>
```

### Frequency Response

![Butterworth frequency response plot](includes/plot_control.svg)

### Comparison

| Design | Passband | Stopband | Notes |
|---|---|---|---|
| Butterworth | Flat | Monotone | General purpose |
| Chebyshev Type I | Ripple ≤ N dB | Steeper | Better rolloff, slight passband distortion |

## Biquad Sections

A `Biquad<T>` implements a single second-order IIR section with transfer function

$$
H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}},
$$

realized in **Direct Form II Transposed** (DFII-T), which is numerically better conditioned than Direct Form I:

```
y[n] = b0*x[n] + s1[n-1]
s1[n] = b1*x[n] - a1*y[n] + s2[n-1]
s2[n] = b2*x[n] - a2*y[n]
```

Coefficients follow the convention `H(z) = (b0 + b1 z⁻¹ + b2 z⁻²) / (1 + a1 z⁻¹ + a2 z⁻²)`.

```rust
use numeris::control::Biquad;

// Construct manually (e.g., from external design tool)
let mut bq = Biquad::new(
    [1.0_f64, 2.0, 1.0],   // [b0, b1, b2]
    [1.0, -1.8, 0.81],     // [a0 (must be 1.0), a1, a2]
);

let y = bq.tick(1.0);
```

`BiquadCascade<T, N>` chains `N` biquad sections in series.

## PID Controller

Discrete-time realization of the PID control law

$$
u(t) = K_p\, e(t) + K_i \int_0^t e(\tau)\,d\tau + K_d\,\frac{d e(t)}{dt},
\qquad e = \text{setpoint} - \text{measurement},
$$

with:

- **Trapezoidal integration** (bilinear, no integrator windup from step changes)
- **Derivative on measurement** (avoids derivative kick on setpoint changes)
- **Optional derivative LPF** (reduces noise amplification)
- **Anti-windup** via back-calculation (clamps integrator when output saturates)
- **Output limits** (hard clamp)

```rust
use numeris::control::Pid;

// Kp=2.0, Ki=0.5, Kd=0.1, sample period dt=0.001 s (1 kHz)
let mut pid = Pid::new(2.0_f64, 0.5, 0.1, 0.001)
    .with_output_limits(-10.0, 10.0)           // clamp output to ±10
    .with_derivative_filter(0.01);              // LPF time constant τ=10ms

let setpoint    = 1.0_f64;
let mut process = 0.0_f64;

for _ in 0..1000 {
    let u = pid.tick(setpoint, process);
    // Simple first-order plant: τ_plant = 0.1 s
    process += 0.001 * (-process + u);
}

assert!((process - setpoint).abs() < 0.01);
```

### API

```rust
// Constructor
fn Pid::new(kp: T, ki: T, kd: T, dt: T) -> Pid<T>

// Builders
fn with_output_limits(self, min: T, max: T) -> Pid<T>
fn with_derivative_filter(self, tau: T)     -> Pid<T>   // τ = filter time constant

// Runtime
fn tick(&mut self, setpoint: T, measurement: T) -> T   // compute output u[k]
fn reset(&mut self)                                      // clear integrator and derivative state
```

### Difference Equations

Let `e[k] = setpoint[k] - measurement[k]`. With trapezoidal integration and derivative on measurement:

```
I[k] = I[k-1] + Ki * dt/2 * (e[k] + e[k-1])     (bilinear integration)
D[k] = -Kd/τ * (measurement[k] - measurement[k-1]) / dt + (1 - dt/τ) * D[k-1]  (filtered derivative)
u[k] = clamp(Kp*e[k] + I[k] + D[k])
```

Anti-windup: if `u` is clamped, the integrator is back-corrected: `I[k] -= Kb * (u_unclamped - u_clamped)`.

## Lead/Lag Compensators

Design continuous-time lead or lag compensators and convert to discrete-time biquad filters via the bilinear transform.

### Lead Compensator

Adds phase lead to improve stability margins. Parameterized by desired maximum phase lead and the frequency where it occurs.

```rust
use numeris::control::{lead_compensator, Biquad};

// 45° phase lead at 50 Hz, unity gain, 1 kHz sample rate
let mut comp = lead_compensator(
    std::f64::consts::FRAC_PI_4, 50.0, 1.0, 1000.0,
).unwrap();
let y = comp.tick(1.0);
```

### Lag Compensator

Boosts low-frequency gain for improved steady-state accuracy without significantly affecting phase margin.

```rust
use numeris::control::{lag_compensator, Biquad};

// 10× DC boost with corner at 5 Hz, 1 kHz sample rate
let mut comp = lag_compensator(10.0, 5.0, 1000.0).unwrap();
let y = comp.tick(1.0);
```

### Bode Plots

![Lead/Lag compensator Bode plots](includes/plot_lead_lag.svg)

## PID Tuning

Automatic PID gain computation from process models or relay experiments.

### FOPDT Model Tuning

Model the process as First-Order Plus Dead Time: `G(s) = K·e^(-Ls) / (τs + 1)`, then apply a tuning rule:

```rust
use numeris::control::{FopdtModel, Pid};

// Identify process: gain=2.5, time constant=3s, dead time=0.5s
let model = FopdtModel::new(2.5, 3.0, 0.5).unwrap();

// Ziegler-Nichols (aggressive, ~25% overshoot)
let zn = model.ziegler_nichols();

// Cohen-Coon (better for large dead time)
let cc = model.cohen_coon();

// SIMC with tau_c = L (balanced performance/robustness)
let simc = model.simc(0.5);

// Apply tuning to a PID controller at 100 Hz
let dt = 0.01;
let mut pid = Pid::new(simc.kp, simc.ki, simc.kd, dt)
    .with_output_limits(-10.0, 10.0);
```

### Ultimate Gain Tuning

From relay or bump-test experiments where you've found the ultimate gain `Ku` and ultimate period `Tu`:

```rust
use numeris::control::ziegler_nichols_ultimate;

let gains = ziegler_nichols_ultimate(10.0_f64, 0.5).unwrap();
// gains.kp = 6.0, gains.ki = 24.0, gains.kd = 0.375
```

### Tuning Rules Summary

| Method | Source | Character |
|---|---|---|
| `ziegler_nichols()` | FOPDT model | Aggressive, ~25% overshoot |
| `cohen_coon()` | FOPDT model | Better for L/τ > 0.5 |
| `simc(tau_c)` | FOPDT model | Adjustable aggressiveness |
| `ziegler_nichols_ultimate()` | Ku, Tu experiment | Classic closed-loop tuning |

## Error Handling

```rust
use numeris::control::ControlError;

match butterworth_lowpass::<f64, 2>(4, 1000.0, 8000.0) {
    Ok(filter) => { /* use filter */ }
    Err(ControlError::InvalidOrder)     => { /* order too large for N sections */ }
    Err(ControlError::InvalidFrequency) => { /* cutoff >= nyquist */ }
    Err(ControlError::InvalidRipple)    => { /* ripple_db <= 0 */ }
}
```
