Skip to content

Tutorial 79: Neuromorphic Control

Spike-domain control theory: PID, Kalman filter, LQR. All controllers use population-coded spike representations. No other SNN library provides control-theory primitives.

Spiking PID Controller

import numpy as np
from sc_neurocore.control import SpikingPID

pid = SpikingPID(Kp=1.0, Ki=0.1, Kd=0.01, n_neurons=10, dt=0.01)

setpoint, measurement = 1.0, 0.0
for step in range(200):
    error = setpoint - measurement
    control = pid.step(error)
    measurement += control * 0.01

# Spike-domain output: population-coded P/I/D channels
rng = np.random.RandomState(42)
spike_output = pid.step_spike(error=0.5, rng=rng)
# shape: (30,) = [P(10), I(10), D(10)]

Spiking Kalman Filter

from sc_neurocore.control import SpikingKalmanFilter

kf = SpikingKalmanFilter(
    n_states=4, n_measurements=2,
    A=np.array([[1,0,0.1,0],[0,1,0,0.1],[0,0,1,0],[0,0,0,1]]),
    H=np.array([[1,0,0,0],[0,1,0,0]]),
)

for t in range(50):
    z = np.array([t*0.1 + np.random.randn()*0.1, t*0.05 + np.random.randn()*0.1])
    state = kf.step(z)

Spiking LQR

from sc_neurocore.control import SpikingLQR

A = np.array([[1.0, 0.1], [0.0, 1.0]])
B = np.array([[0.0], [0.1]])
lqr = SpikingLQR(A=A, B=B)

x = np.array([1.0, 0.0])
for t in range(100):
    u = lqr.control(x)
    x = A @ x + B @ u
# x converges to origin

API Reference

sc_neurocore.control.controllers

Spike-domain control: PID, Kalman filter, LQR.

All controllers use population-coded spike representations. Gains are synaptic weights, integration is membrane dynamics. No SNN library provides control-theory primitives.

Stagsted et al. 2020 (RSS) — spiking PID on Loihi

SNN-LQR-EMSIF (Nature Scientific Reports 2025)

SpikingPID

Population-coded PID controller.

Error → rate-coded spike populations → P/I/D populations → output current. Gains are synaptic weights.

Parameters

Kp, Ki, Kd : float PID gains (encoded as synaptic weights). n_neurons : int Population size per channel. dt : float Timestep.

Source code in src/sc_neurocore/control/controllers.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
class SpikingPID:
    """Population-coded PID controller.

    Error → rate-coded spike populations → P/I/D populations →
    output current. Gains are synaptic weights.

    Parameters
    ----------
    Kp, Ki, Kd : float
        PID gains (encoded as synaptic weights).
    n_neurons : int
        Population size per channel.
    dt : float
        Timestep.
    """

    def __init__(
        self,
        Kp: float = 1.0,
        Ki: float = 0.1,
        Kd: float = 0.01,
        n_neurons: int = 10,
        dt: float = 0.01,
    ):
        self.Kp = Kp
        self.Ki = Ki
        self.Kd = Kd
        self.n_neurons = n_neurons
        self.dt = dt
        self._integral = 0.0
        self._prev_error = 0.0

    def step(self, error: float) -> float:
        """Compute PID output for one timestep.

        Parameters
        ----------
        error : float
            Setpoint - measurement.

        Returns
        -------
        float — control output
        """
        self._integral += error * self.dt
        derivative = (error - self._prev_error) / self.dt if self.dt > 0 else 0.0
        self._prev_error = error
        return self.Kp * error + self.Ki * self._integral + self.Kd * derivative

    def step_spike(self, error: float, rng: np.random.RandomState | None = None) -> np.ndarray:
        """Compute PID output as spike population.

        Returns binary spike vector of shape (3 * n_neurons,) representing
        [P_population, I_population, D_population].
        """
        if rng is None:
            rng = np.random.RandomState(0)
        output = self.step(error)

        # Population-code each component
        p_rate = np.clip(abs(self.Kp * error) / 10, 0, 1)
        i_rate = np.clip(abs(self.Ki * self._integral) / 10, 0, 1)
        d_rate = np.clip(abs(self.Kd * (error - self._prev_error)) / 10, 0, 1)

        p_spikes = (rng.random(self.n_neurons) < p_rate).astype(np.int8)
        i_spikes = (rng.random(self.n_neurons) < i_rate).astype(np.int8)
        d_spikes = (rng.random(self.n_neurons) < d_rate).astype(np.int8)

        return np.concatenate([p_spikes, i_spikes, d_spikes])

    def reset(self):
        self._integral = 0.0
        self._prev_error = 0.0

step(error)

Compute PID output for one timestep.

Parameters

error : float Setpoint - measurement.

Returns

float — control output

Source code in src/sc_neurocore/control/controllers.py 
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def step(self, error: float) -> float:
    """Compute PID output for one timestep.

    Parameters
    ----------
    error : float
        Setpoint - measurement.

    Returns
    -------
    float — control output
    """
    self._integral += error * self.dt
    derivative = (error - self._prev_error) / self.dt if self.dt > 0 else 0.0
    self._prev_error = error
    return self.Kp * error + self.Ki * self._integral + self.Kd * derivative

step_spike(error, rng=None)

Compute PID output as spike population.

Returns binary spike vector of shape (3 * n_neurons,) representing [P_population, I_population, D_population].

Source code in src/sc_neurocore/control/controllers.py 
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def step_spike(self, error: float, rng: np.random.RandomState | None = None) -> np.ndarray:
    """Compute PID output as spike population.

    Returns binary spike vector of shape (3 * n_neurons,) representing
    [P_population, I_population, D_population].
    """
    if rng is None:
        rng = np.random.RandomState(0)
    output = self.step(error)

    # Population-code each component
    p_rate = np.clip(abs(self.Kp * error) / 10, 0, 1)
    i_rate = np.clip(abs(self.Ki * self._integral) / 10, 0, 1)
    d_rate = np.clip(abs(self.Kd * (error - self._prev_error)) / 10, 0, 1)

    p_spikes = (rng.random(self.n_neurons) < p_rate).astype(np.int8)
    i_spikes = (rng.random(self.n_neurons) < i_rate).astype(np.int8)
    d_spikes = (rng.random(self.n_neurons) < d_rate).astype(np.int8)

    return np.concatenate([p_spikes, i_spikes, d_spikes])

SpikingKalmanFilter

Spike-domain Kalman filter for state estimation.

State prediction and correction using LIF-based integration. Kalman gain encoded as synaptic weight matrix.

Parameters

n_states : int State dimension. n_measurements : int Measurement dimension. A : ndarray State transition matrix. H : ndarray Observation matrix. Q : ndarray Process noise covariance. R : ndarray Measurement noise covariance.

Source code in src/sc_neurocore/control/controllers.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
class SpikingKalmanFilter:
    """Spike-domain Kalman filter for state estimation.

    State prediction and correction using LIF-based integration.
    Kalman gain encoded as synaptic weight matrix.

    Parameters
    ----------
    n_states : int
        State dimension.
    n_measurements : int
        Measurement dimension.
    A : ndarray
        State transition matrix.
    H : ndarray
        Observation matrix.
    Q : ndarray
        Process noise covariance.
    R : ndarray
        Measurement noise covariance.
    """

    def __init__(
        self,
        n_states: int,
        n_measurements: int,
        A: np.ndarray | None = None,
        H: np.ndarray | None = None,
        Q: np.ndarray | None = None,
        R: np.ndarray | None = None,
    ):
        self.n_states = n_states
        self.n_measurements = n_measurements
        self.A = A if A is not None else np.eye(n_states)
        self.H = H if H is not None else np.eye(n_measurements, n_states)
        self.Q = Q if Q is not None else np.eye(n_states) * 0.01
        self.R = R if R is not None else np.eye(n_measurements) * 0.1
        self.x = np.zeros(n_states)
        self.P = np.eye(n_states)

    def predict(self) -> np.ndarray:
        """Predict step: x = A @ x, P = A @ P @ A^T + Q."""
        self.x = self.A @ self.x
        self.P = self.A @ self.P @ self.A.T + self.Q
        return self.x.copy()

    def update(self, z: np.ndarray) -> np.ndarray:
        """Update step with measurement z."""
        S = self.H @ self.P @ self.H.T + self.R
        K = self.P @ self.H.T @ np.linalg.inv(S)
        innovation = z - self.H @ self.x
        self.x = self.x + K @ innovation
        self.P = (np.eye(self.n_states) - K @ self.H) @ self.P
        return self.x.copy()

    def step(self, z: np.ndarray) -> np.ndarray:
        """Predict + update in one call."""
        self.predict()
        return self.update(z)

    def reset(self):
        self.x = np.zeros(self.n_states)
        self.P = np.eye(self.n_states)

predict()

Predict step: x = A @ x, P = A @ P @ A^T + Q.

Source code in src/sc_neurocore/control/controllers.py 
139
140
141
142
143
def predict(self) -> np.ndarray:
    """Predict step: x = A @ x, P = A @ P @ A^T + Q."""
    self.x = self.A @ self.x
    self.P = self.A @ self.P @ self.A.T + self.Q
    return self.x.copy()

update(z)

Update step with measurement z.

Source code in src/sc_neurocore/control/controllers.py 
145
146
147
148
149
150
151
152
def update(self, z: np.ndarray) -> np.ndarray:
    """Update step with measurement z."""
    S = self.H @ self.P @ self.H.T + self.R
    K = self.P @ self.H.T @ np.linalg.inv(S)
    innovation = z - self.H @ self.x
    self.x = self.x + K @ innovation
    self.P = (np.eye(self.n_states) - K @ self.H) @ self.P
    return self.x.copy()

step(z)

Predict + update in one call.

Source code in src/sc_neurocore/control/controllers.py 
154
155
156
157
def step(self, z: np.ndarray) -> np.ndarray:
    """Predict + update in one call."""
    self.predict()
    return self.update(z)

SpikingLQR

Spike-domain Linear Quadratic Regulator.

Computes optimal gain K from system matrices (A, B, Q, R). Control law: u = -K @ x. Weights derived analytically.

Parameters

A : ndarray (n, n) — state transition B : ndarray (n, m) — control input Q : ndarray (n, n) — state cost R : ndarray (m, m) — control cost

Source code in src/sc_neurocore/control/controllers.py 
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
class SpikingLQR:
    """Spike-domain Linear Quadratic Regulator.

    Computes optimal gain K from system matrices (A, B, Q, R).
    Control law: u = -K @ x. Weights derived analytically.

    Parameters
    ----------
    A : ndarray (n, n) — state transition
    B : ndarray (n, m) — control input
    Q : ndarray (n, n) — state cost
    R : ndarray (m, m) — control cost
    """

    def __init__(
        self, A: np.ndarray, B: np.ndarray, Q: np.ndarray | None = None, R: np.ndarray | None = None
    ):
        n = A.shape[0]
        m = B.shape[1]
        self.A = A
        self.B = B
        self.Q = Q if Q is not None else np.eye(n)
        self.R = R if R is not None else np.eye(m)
        self.K = self._solve_dare()

    def _solve_dare(self, max_iter: int = 200) -> np.ndarray:
        """Solve discrete algebraic Riccati equation iteratively."""
        P = self.Q.copy()
        for _ in range(max_iter):
            K = np.linalg.solve(
                self.R + self.B.T @ P @ self.B,
                self.B.T @ P @ self.A,
            )
            P_new = self.Q + self.A.T @ P @ (self.A - self.B @ K)
            if np.allclose(P, P_new, atol=1e-10):
                break
            P = P_new
        return np.linalg.solve(
            self.R + self.B.T @ P @ self.B,
            self.B.T @ P @ self.A,
        )

    def control(self, x: np.ndarray) -> np.ndarray:
        """Compute optimal control: u = -K @ x."""
        return -self.K @ x

    @property
    def gain_matrix(self) -> np.ndarray:
        return self.K.copy()

control(x)

Compute optimal control: u = -K @ x.

Source code in src/sc_neurocore/control/controllers.py 
206
207
208
def control(self, x: np.ndarray) -> np.ndarray:
    """Compute optimal control: u = -K @ x."""
    return -self.K @ x