Hello CVX community,
I am implementing the data-driven LQR formulation from De Persis & Tesi (2020) “Formulas for Data-Driven Control: Stabilization, Optimality, and Robustness” (IEEE TAC) and am getting the following warning from SeDuMi, after which the solve returns “Failed”.
“The coefficient matrix is not full row rank, numerical problems may occur. SeDuMi 1.3.7 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.”
Below is a minimum reproducible example snippet of code.
clear
clc
rng(42) % Fixed seed, disable to reproduce stoachastic failure
%% System Definition
A = [-0.7833, 0.0668, 0.0004, 0;
0.3752, -0.0668, 0, 0;
0.0196, 0, -0.0004, 0;
0.4560, 0, 0, -0.4560];
B = [0.2342; 0; 0; 0];
C50 = 6.33;
gamma = 2.24;
E_0 = 98.8;
E_max = 94.1;
% Additional system to test with
% A = [-0.7943, 0.0665, 0.0004, 0;
% 0.3977, -0.0665, 0, 0;
% 0.0196, 0, -0.0004, 0;
% 0.4560, 0, 0, -0.4560];
% B = [0.2342; 0; 0; 0];
%C50 = 6.76;
%gamma = 4.29;
%E_0 = 98.6;
%E_max = 86.0;
% Find equilibrium point and collect data about equilibrium
[x_eq, u_eq] = calc_eq(A, B, 50, C50, gamma, E_0, E_max);
%% Persistently Exciting Data Collection
% Minimum samples required for PE: N_data >= m + n = 5
dt_sim = 0.01; % Time step for simulation
dt_data = 0.5; % Data collection rate
N_data = 6; % Minimum number (+1 for equilibrium states)
f_ode = @(x,u) A*x + B*u;
x_data(:,1) = x_eq;
x_current = x_data(:,1);
u_data(:,1) = u_eq;
delta_u = 2;
for l = 1:N_data
du = delta_u * randn(); % Excitation signal
u_data(l+1) = u_eq + du;
for i = 1:dt_data/dt_sim
x_current = rk4(f_ode, x_current, u_data(l), dt_sim);
end
x_data(:,l+1) = x_current;
end
%% Solve for data driven LQR
% Theory: De Persis & Tesi (2020) "Formulas for Data-Driven Control"
Q = diag([10,1,1,500]);
r = 1;
K_dlqr = findK_data(x_data(:,2:end),u_data(2:end),Q,r);
%% Functions %%
function [x_eq, u_eq] = calc_eq(A, B, BIS_target, C50, gamma, E0, Emax)
% Calculate equilibrium effect-site concentration
Ce_eq = C50 * ((E0 - BIS_target) / (Emax - E0 + BIS_target))^(1/gamma);
% Calculate steady-state gain from u to x
tmp = -A\B;
u_eq = Ce_eq/tmp(4); % Calculate equilibrium infusion rate
x_eq = tmp*u_eq; % Calculate full equilibrium state
end
function K_lqr = findK_data(X,U,Qx,R)
n = size(X,1);
m = size(U,1);
T = size(X,2) - 1;
L = 1;
R_sqrt = sqrtm(R); % Tuning parameters
X0 = X(:,1:end-1);
X1 = X(:,2:end);
U0 = U(:,1:end-1);
% Enforce full rank condition
UX = cat(1,U0,X0);
cond(UX)
if rank(UX) ~= (m+n)
error("Signal is NOT persistently exciting")
end
% Solve for LQR optimal gain. NOTE: Numeric instability is quite often seen here
cvx_clear
cvx_begin sdp
cvx_precision default
variable Q_lqr(T/L, n)
variable Y(m, m) symmetric
variable P_lqr(n, n) symmetric
minimize(trace(Qx*X0*Q_lqr) + trace(Y))
subject to
% Schur complement for cost
[Y, R_sqrt*U0*Q_lqr;
Q_lqr'*U0'*R_sqrt', P_lqr] >= 0;
% Schur complement for Lyapunov stability
[P_lqr - eye(n), X1*Q_lqr;
Q_lqr'*X1', P_lqr] >= 0;
P_lqr == X0*Q_lqr;
cvx_end
% Check CVX status
if ~strcmp(cvx_status, 'Solved')
warning('CVX optimization failed with status: %s', cvx_status);
end
K_lqr = U0 * Q_lqr / (X0*Q_lqr); % Recover controller gain
end
function x_next = rk4(f,x,u,dt)
% Fourth-order Runge-Kutta
k1 = f(x,u);
k2 = f((x + 0.5*dt*k1), u);
k3 = f((x + 0.5*dt*k2), u);
k4 = f((x + dt*k3), u);
x_next = x + (dt/6) * (k1 + 2*k2 + 2*k3 + k4);
end
I believe the implementation follows the formulation correctly, but I acknowledge there may still be subtle issues in how the SDP is posed. This is one possible cause I have been considering:
- The system is stiff (eigenvalues span
-0.0004to-0.81), givingcond([U0; X0]) ≈ 1e7-1e10even when full rank. Collecting more data would help numerically, but this formulation is appealing precisely because it operates at the theoretical minimum data length, so I would prefer not to rely on that as the fix. Furthermore, my current sampling time is relatively high, which is costing controller performance.
I am particularly interested in whether this issue is inherent to the SDP formulation itself. More generally, are there theoretical or practical considerations (e.g., redundancy in constraints, scaling issues, or solver limitations) that could explain why SeDuMi reports a lack of full row rank?
Thank you for your time and patience.