SeDuMi "coefficient matrix not full row rank" for data-driven LQR SDP

1

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:

  1. The system is stiff (eigenvalues span -0.0004 to -0.81), giving cond([U0; X0]) ≈ 1e7-1e10 even 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.