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

Question

Maureen 2024-3-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2095551-why-do-i-receive-this-error-error-using-odearguments-odes_lin-must-return-a-column-vector

评论： Sam Chak 2024-3-18

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

7 个评论
显示 5更早的评论隐藏 5更早的评论

Maureen 2024-3-17

在 MATLAB Online 中打开

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

Torsten 2024-3-17

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 2024-3-18

编辑：Stephen23 2024-3-18

在 MATLAB Online 中打开

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 2024-3-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2095551-why-do-i-receive-this-error-error-using-odearguments-odes_lin-must-return-a-column-vector#answer_1426686

移动：Cris LaPierre 2024-3-18

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;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Sam Chak 2024-3-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2095551-why-do-i-receive-this-error-error-using-odearguments-odes_lin-must-return-a-column-vector#answer_1426811

在 MATLAB Online 中打开

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Sam Chak 2024-3-18

在 MATLAB Online 中打开

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.