Why do I Receive this error: Error using odearguments ODES_LIN must return a column vector.

I created an H-infinity controller for a system and am attempting implement the controller in a numerical simulation to plot the system response and control output, but receive the following error: Error using odearguments ODES_LIN must return a column vector.

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in proposalGPBode (line 149)

[t,x_out] = ode45(@ODEs_lin,t_span,IC-[x_bar;0*xc0],options,const,Ap,Bp,Cp,Ac,Bc,Cc);

I've included Ap,Bp,Cp, and Dp in the code below, but just have the matrix sizes of Ac, Bc, and Dc here.

Ac: 7x7 double

Bc:7x1 double

Cc:1x7 double

Any ideas on why I am not returning a column vector?

u_bar = -b*sqrt(-1/(c*a));
x_bar = [0;b/c];

Ap=[0 b;0 -c];
Bp=[2*a*u_bar; 0];
Cp=[1 0];
Dp=0;

%% Initial conditions
xc0 = zeros(size(Ac,1),1);
IC = [0; 0.1; xc0]; % Initial state [x1; x2]

% Simulation time.
t0 = 0; % initial time (s)
t_max = 5; % final time (s)
t_div = 1001; % number of steps to divide the time series into when plotting
t_span = linspace(t0,t_max,t_div); % define simulation time

% Simulation options.
 options = odeset('AbsTol',1e-5,'RelTol',1e-5); % This changes the 
 % integration tolerence.  Smaller values of AbsTol and RelTol (e.g., 1e-12
 % versus 1e-6) will give more accurate numerical integration solutions.

 % Numerical solve ODE

    [t,x_out] = ode45(@ODEs_lin,t_span,IC-[x_bar;0*xc0],options,const,Ap,Bp,Cp,Ac,Bc,Cc);


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Linear Differential Equation to Numerically Integrate

function  [dot_x] = ODEs_lin(t,x,const,Ap,Bp,Cp,Ac,Bc,Cc)

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate Control Input

delta_x1 = x(1);
delta_x2 = x(2);
delta_x=[delta_x1;delta_x2];
xc = x(3:2+size(Ac,1));

[d,n] = ExogenousSignals(t,const);

 y = -Cp*delta_x - n;
delta_u=Cc*xc;

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dynamics

dot_delta_x = Ap*delta_x + Bp*d+ Bp*delta_u;
dot_xc = Ac*xc + Bc*y;
dot_x = [dot_delta_x;dot_xc];
end

function  [d,n] = ExogenousSignals(t,const)

d = sum(const.dist_amp.*sin(const.dist_freq*2*pi*t));
n = sum(const.noise_amp.*sin(const.noise_freq*2*pi*t));

end

[Torsten]

Please insert the line

size(dot_x)

at the end of function "ODEs_lin" and tell us what MATLAB prints out.

[Maureen]

size(dot_x)= 9 5

[Torsten]

dot_x must be a vector of size (nx x 1) where nx is the number of equations you are trying to solve. So in your case maybe 45 x 1.

We cannot run your code to test it because variables are not defined.

[Maureen]

clc
clear all
close all

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options

Nonlin_bool = 0; % Choose 0 for linear simulation and 1 for nonlinear simulation;
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constants 
a=-1; % Example value for the parameter alpha
c = 1; % Example value for the parameter c
b=1; % Example value for the parameter beta

% Define noise properties for simulations
const.dist_freq = 5*rand(2,5);
const.dist_amp = rand(2,5)/5;
const.noise_freq = 50*randn(3,5)+100;
const.noise_amp = diag([0.01;0.01;1*pi/180])*rand(3,5)/5;


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Controller Design

% Linearization
u_bar = -b*sqrt(-1/(c*a));
x_bar = [0;b/c];

Ap=[0 b;0 -c];% eqn 1 & 2:B
Bp=[2*a*u_bar; 0];
Cp=[1 0];
Dp=0;
Gp = ss(Ap,Bp,Cp,Dp);

% Define weights

s = tf('s');

% Wd is 1st-order low-pass filter
d_cutoff = 2*pi*5; %5 Hz cut off
d_mag_max = 1e-6; % Maximum of 1 uN disturbance expected
Wd = d_mag_max*d_cutoff/(s+d_cutoff);

% Wn is 1st-order high-pass filter with feedthough
n_low_cutoff = 2*pi*30000;
n_high_cutoff = 2*pi*1000000;
n_mag_pos_max = 1.4e-3; % Maximum of 1.4 mV noise expected pmt 1 & pmt 2
n_deg = 1; % order/degree of high-pass filter
Wn = n_high_cutoff/n_low_cutoff*(s+n_low_cutoff)/(s+n_high_cutoff)*n_mag_pos_max;

% Wr is 1st-order high-pass filter with feedthough
r_low_cutoff = 2*pi*100;
r_high_cutoff = 2*pi*10000;
r_mag_pos_max = 0.01e-3; % Maximum of 1.4 mV noise expected pmt 1 & pmt 2
r_deg = 1; % order/degree of high-pass filter
Wr = r_high_cutoff/r_low_cutoff*(s+r_low_cutoff)/(s+r_high_cutoff)*r_mag_pos_max;

% We is 1st-order low-pass filter with feedthrough
e_low_cutoff = 2*pi*100;
e_high_cutoff = 2*pi*10000;
e_mag_pos_max = 0.01; % Maximum of 0.01 V error at low freq.
e_deg = 1; % order/degree of low-pass filter
We = e_low_cutoff^e_deg/e_high_cutoff^e_deg*(s+e_high_cutoff)^e_deg/(s+e_low_cutoff)^e_deg*(1/e_mag_pos_max);

We_const = ss([],[],[],1/e_mag_pos_max);

% Wu is 1st-order high-pass filter with feedthrough
u_low_cutoff = 2*pi*1;
u_high_cutoff = 8000;
u_mag_max = 100e-3; % Maximum of 100 mN input at low freq.
u_deg = 2; % order/degree of high-pass filter
Wu = u_high_cutoff^u_deg/u_low_cutoff^u_deg/u_mag_max*(s+u_low_cutoff)^u_deg/(s+u_high_cutoff)^u_deg;

Wu_const = ss([],[],[],1/u_mag_max);

% Define weighted generalized plant

Gp.u = 'u_p';
Gp.y = 'y_p';

Wd.u = 'd';
Wd.y = 'd_weight';
Wn.u = 'n';
Wn.y = 'n_weight';
Wr.u = 'r';
Wr.y = 'r_weight';
We.u = 'e';
We.y = 'e_weight';
Wu.u = 'u';
Wu.y = 'u_weight';

We_const.u = 'e';
We_const.y = 'e_const';
Wu_const.u = 'u';
Wu_const.y = 'u_const';

input = {'d', 'n','r', 'u'};
output = {'e_const', 'u_const', 'y'};
output_weighted = {'e_weight', 'u_weight', 'y'};


Sum1 = sumblk('u_p = d_weight + u'); 
Sum2 = sumblk('y = r_weight- y_p - n_weight'); 
Sum3 = sumblk('e = r-y_p');

P_gen_weight_Hinf = connect(Gp, Wd, Wn, Wr,We, Wu, Sum1, Sum2, Sum3, input, output_weighted); 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% H2-Optimal and Hinf-Optimal Controller Design

[K_Hinf,CL_Hinf,GAM_Hinf] = hinfsyn(P_gen_weight_Hinf,1,1);

K_sim = K_Hinf;

% Extract state space matrices from the controller for the simulation
Ac = K_sim.a;
Bc = K_sim.b;
Cc = K_sim.c;
Dc = K_sim.d;

% Example of Bode plots
figure
bode(K_Hinf)

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Initial conditions
xc0 = zeros(size(Ac,1),1);
IC = [0; 0.1; xc0]; % Initial state [x1; x2]

% Simulation time.
t0 = 0; % initial time (s)
t_max = 5; % final time (s)
t_div = 1001; % number of steps to divide the time series into when plotting
t_span = linspace(t0,t_max,t_div); % define simulation time

% Simulation options.
 options = odeset('AbsTol',1e-5,'RelTol',1e-5); % This changes the 
 % integration tolerence.  Smaller values of AbsTol and RelTol (e.g., 1e-12
 % versus 1e-6) will give more accurate numerical integration solutions.

 % Numerical solve ODE
 if Nonlin_bool == 1
    [t,x_out] = ode45(@ODEs,t_span,IC,options,const,x_bar,u_bar,Cp,Ac,Bc,Cc);
 else
    [t,x_out] = ode45(@ODEs_lin,t_span,IC-[x_bar;0*xc0],options,const,Ap,Bp,Cp,Ac,Bc,Cc);
 end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Post Processing

for lv1 = 1:length(x_out)
    % Note, I like to use ``lv1'' for my loop counter, which stands for
    % ``loop-varibale 1'' because the letters i and j are reserved as
    % imaginary numbers.
        
    if Nonlin_bool == 1
        delta_x(:,lv1)  = x_out(lv1,1:2)' - x_bar;  
    else
        delta_x(:,lv1)  = x_out(lv1,1:2)';
    end
    
    [d,n] = ExogenousSignals(t(lv1),const);
    
    x(:,lv1) = delta_x(:,lv1) + x_bar;
    xc(:,lv1) =x_out(lv1,3:2+size(Ac,1))';
    y = -Cp*delta_x(:,lv1) - n;
        
    delta_u(:,lv1)  = Cc*xc(:,lv1);
    u(:,lv1)        = u_bar + delta_u(:,lv1);
    
end


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot data

% Font size, line size, and line width. 
font_size = 15;
line_size = 15;
line_width = 2;

% Notice how the axes for each plot are labeled with the appropriate units

% System response
figure
subplot(2,1,1)
plot(t,x(1,:),'Linewidth',line_width);
hold on
title('System Trajectories vs Time')
xlabel('Time (us)','fontsize',font_size,'Interpreter','latex');
ylabel('X1 (um)','fontsize',font_size,'Interpreter','latex');
set(gca,'XMinorGrid','off','GridLineStyle','-','FontSize',line_size)
grid on
subplot(2,1,2)
plot(t,x(2,:)*180/pi,'Linewidth',line_width);
hold on
xlabel('Time (us)','fontsize',font_size,'Interpreter','latex');
ylabel('X2','fontsize',font_size,'Interpreter','latex');
set(gca,'XMinorGrid','off','GridLineStyle','-','FontSize',line_size)
grid on


% % Control Inputs
figure
plot(t,u(1,:),'Linewidth',line_width);
hold on
title('Control Input vs Time');
xlabel('Time (us)','fontsize',font_size,'Interpreter','latex');
ylabel('u(kV)','fontsize',font_size,'Interpreter','latex');
set(gca,'XMinorGrid','off','GridLineStyle','-','FontSize',line_size)
grid on

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Linear Differential Equation to Numerically Integrate

function  [dot_x] = ODEs_lin(t,x,const,Ap,Bp,Cp,Ac,Bc,Cc)

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate Control Input

delta_x1 = x(1);
delta_x2 = x(2);
delta_x=[delta_x1;delta_x2];
xc = x(3:2+size(Ac,1));

[d,n] = ExogenousSignals(t,const);

 y = -Cp*delta_x - n;
delta_u=Cc*xc;

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dynamics

dot_delta_x = Ap*delta_x + Bp*d+ Bp*delta_u;
dot_xc = Ac*xc + Bc*y;
dot_x = [dot_delta_x;dot_xc];

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nonlinear Differential Equation to Numerically Integrate

function  [dot_x] = ODEs(t,x,const,x_bar,u_bar,Cp,Ac,Bc,Cc)
% 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Calculate Control Input

delta_x = x(1:2)-x_bar;
xc = x(3:2+size(Ac,1));

[dist,n] = ExogenousSignals(t,const)

y = -Cp*delta_x - n;
delta_u = Cc*xc;
u = u_bar + delta_u;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dynamics
a=-1; % Example value for the parameter alpha
c = 1; % Example value for the parameter c
b=1; % Example value for the parameter beta

% First, extract states in a convenient form. 
x1 = x(1);
x2 = x(2);

% Form dot_x = f(x,u) system.
x1_dot=x2*b+a*u*u;
x2_dot=-c*x2 + b;
dot_x_plant = [x1_dot;x2_dot];
dot_xc = Ac*xc + Bc*y;
dot_x = [dot_x_plant;dot_xc];
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  [d,n] = ExogenousSignals(t,const)

d = sum(const.dist_amp.*sin(const.dist_freq*2*pi*t));
n = sum(const.noise_amp.*sin(const.noise_freq*2*pi*t));

end

[Maureen]

Here is the code: any ideas on how to reformat dot_x to make it a column vector?

[Torsten]

Your vector of initial conditions IC-[x_bar;0*xc0] has size 9x1. Thus dot_x must be a vector of size 9x1. But you produce a matrix of size 9x5. I don't know what's going wrong in your computations to get dot_x.

[Stephen23]

Note that you are passing extra arguments using an approach that has not been documented for twenty years. Your approach has been deprecated ever since anonymous functions were introduced in R14 exactly twenty years ago.

I strongly recommend that you parameterize your functions using the documented approaches:

https://www.mathworks.com/help/matlab/math/parameterizing-functions.html

For example, using an anonymous function, something like this:

fnh = @(t,x) ODEs(t,x, const,x_bar,u_bar,Cp,Ac,Bc,Cc);
[t,x_out] = ode45(fnh, t_span, IC, options);

Answer 1

[Maureen]

I think I figured it out. When I redefined these noise and disturbance equations, it runs. Thank you for your help!

const.dist_freq = 5*rand(2,5);

const.dist_amp = rand(2,5)/5;

const.noise_freq = 50*randn(3,5)+100;

const.noise_amp = diag([0.01;0.01;1*pi/180])*rand(3,5)/5;

Answer 2

[Sam Chak]

Hi @Maureen

It seems that you intend to generate the exogenous disturbance and noise signals as the sum of sine waves. Specifically, the disturbance signal comprises the sum of two sine waves, while the noise signal consists of the sum of three sine waves. Thus, only one noisy vector 'n' affects the output, and the matrix 'Bp' in the process plant Gp can only accommodate a single disturbance vector 'd'.

%% Settings for Exogenous Signals
num_d            = 2;   % number of disturbing sine waves
num_n            = 3;   % number of noisy sine waves
const.dist_freq  = 5*rand(num_d,1);
const.dist_amp   = rand(num_d,1)/5;
const.noise_freq = 50*randn(num_n,1)+100;
const.noise_amp  = diag([0.01;0.01;1*pi/180])*rand(num_n,1)/5;

%% Test
t       = 0:0.1:5;
[d, n]  = ExogenousSignals(t, const)

d = 1×51

0 0.0929 -0.1014 0.0339 0.0498 -0.1012 0.0815 0.0115 -0.1009 0.1000 -0.0188 -0.0655 0.1002 -0.0673 -0.0235 0.1054 -0.0965 0.0049 0.0802 -0.0986 0.0511 0.0363 -0.1068 0.0901 0.0078 -0.0929 0.0967 -0.0339 -0.0500 0.1056

n = 1×51

0 0.0000 -0.0022 -0.0007 -0.0027 -0.0026 -0.0020 -0.0042 -0.0014 -0.0041 -0.0019 -0.0023 -0.0023 -0.0006 -0.0012 -0.0002 0.0012 -0.0003 0.0032 0.0005 0.0035 0.0024 0.0025 0.0037 0.0021 0.0030 0.0026 0.0011 0.0027 -0.0004

%% Function for creating Exogenous disturbance and noise
function  [d, n] = ExogenousSignals(t, const)
    d = sum(const.dist_amp.*sin(const.dist_freq*2*pi*t));
    n = sum(const.noise_amp.*sin(const.noise_freq*2*pi*t));
end

[Sam Chak]

The synthesized 7th-order H-infinity () controller will generate the following trajectories for the 1st-order Integrator Plant.

%% Constants 
a   =-1; % Example value for the parameter alpha
b   = 1; % Example value for the parameter beta
c   = 1; % Example value for the parameter c

%% Process Plant
u_bar   = -b*sqrt(-1/(c*a));
x_bar   = [0;b/c];
Ap      = [0  b; 
           0 -c];   % eqn 1 & 2:B
Bp      = [2*a*u_bar; 
           0];
Cp      = [1, 0];
Dp      =  0;
Gp      = ss(Ap, Bp, Cp, Dp);
Gp      = tf(Gp)

Gp =
 
  2
  -
  s
 
Continuous-time transfer function.