Error: Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Hello, I'm fairly new to MATLAB and I am unsure as to why my code came up with this error. This was part of the error message:

Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Error in evalcontrol (line 9)
u(:,ii)=-sys.r\sys.b'*reshape(deval(solp,Tf-t),length(x),length(x))*x;

Error in ACS420_3 (line 37)
[u,J,Jf]=evalcontrol(tspan,sys_name,solx,solp,Tf);

My parameters are as follows, and are in a function called imresp.m:

function sys=imresp(t,x)
a11=1; a12=1; a23=1; a31=1; a42=1; b2=1; b3=1;
b1=-1; b4=-1;
a22=3;
a32=1.5;
a33=0.5; a41=0.5;
f11=1; f44=1; q11=1; q44=1; r11=1; r22=1; 
x2bar=2; x3bar=a31*x2bar/a32;
if nargin==0; sys.f = [f11 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 f44]; sys.xe = [0 x2bar x3bar 0]';
else
sys.a = [((a11-a12*x3bar)-a12*x(3)) 0 0 0;...
    a21(x(4))*a22*(x(3)+x3bar) -a23 0 0;...
    -a33*x3bar a31 -(a32+a33*x(1)) 0;...
    a41 0 0 -a42];
sys.b = [b1 0 0 0; 0 b2 0 0]';
sys.q = [q11 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 q44];
sys.r = [r11 0; 0 r22];
end

and evalcontrol.m is:

function [u,J,Jf]=evalcontrol(tspan,sys_name,solx,solp,Tf)
xf=deval(solx,Tf);
tmp=zeros(1,length(tspan));
u=zeros(1,length(tspan));
for ii=1:length(tspan)
    t=tspan(ii);
    x=deval(solx,t);
    sys=feval(sys_name,t,x);
    u(:,ii)=-sys.r\sys.b'*reshape(deval(solp,Tf-t),length(x),length(x))*x;
    tmp(ii)=x'*sys.q*x+u(:,ii)'*sys.r*u(:,ii);
end
J=cumtrapz(tspan,tmp)/2; % cumtrapz - cumulative trapezoidal numerical integration
sys=imresp2;
Jf=0.5*xf'*sys.f*xf;

Help is much appreciated. Thank you!

[Walter Roberson]

It is not clear how your are invoking evalcontrol ?

The initialization code you give does not have a definition for x.

[chromaclouds]

sorry for the lack of clarity - evalcontrol.m is invoked in a main script called ACS420_3.m.

IC = [1.5 2 2.57 3];
sys_name='imresp';

imresp.m has all the parameters I stated above (values of aij, bij, fij, qij, rij, x2bar, x3bar, sys.f, sys.xe, sys.a, sys.b, sys.q, sys.r).

---

also in the main script ACS420_3.m:

xe=sys.xe;
xo=[0 0 0 0]';
solx=[];solp=[];
xo(1)=IC(ii);solx=xo; % First iteration uses fixed value of xo.
options=[];
dT = 0.01;
To = 0; Tf = 10;
tspan = To:dT:Tf;

in ACS420_3.m, while the termination value has not been reached, the following code runs:

solp=ode23s(rhandle,[0 Tf],po,options,sys_name,solx,Tf); % solve riccati ode
solx=ode23s(dhandle,[0 Tf],xo,options,sys_name,solx,solp,Tf); % solve dynamics ode
[u,J,Jf]=evalcontrol(tspan,sys_name,solx,solp,Tf); % compute control signal

where

rhandle=@riccati;

and riccati.m is a function:

function dp=riccati(t,p,sys_name,solx,Tf)
n=sqrt(length(p));
P=reshape(p,n,n);
if ~isstruct(solx);
    sys=feval(sys_name,t,solx); % first iteration sets to fixed value for all t, i.e. the IC.
else
    sys=feval(sys_name,t,deval(solx,Tf-t));
end
dP=sys.q + P*sys.a + sys.a'*P - P*sys.b/sys.r*sys.b'*P;
dp=dP(:);

and solx is where dhandle is:

dhandle=@dynamics;

and dynamics.m is a function:

function dx=dynamics(t,x,sys_name,solx,solp,Tf)
if ~isstruct(solx);
    sys=feval(sys_name,t,solx); % first iteration sets to fixed value for all t, i.e. the IC.
else
    sys=feval(sys_name,t,deval(solx,Tf-t));
end
k=-sys.r\sys.b'*reshape(deval(solp,Tf-t),length(x),length(x));
dx=(sys.a+sys.b*k)*x;

[Walter Roberson]

I think we need all of ACS420_3.m .

You posted riccati and dynamics and imresp and evalcontrol already; if there is anything else beyond those and ACS420_3 then you should post the others as well.

[chromaclouds]

alright, ACS420_3.m is below.

% runimmune CL  solves closed-loop non-linear optimal control for pathogen attack  
% housekeeping
clear; close all
%initial concentration of pathogen
IC = [1.5 2 2.57 3];
sys_name='imresp2'; % define system name
sys=imresp2; % get necessary system info
po=sys.f(:); % vectorize the F matrix
xe=sys.xe;   % this is the equilibrium value so we can re-construct actual concs.
options=[];  % options for ode solver (leaves at default)
dhandle=@dynamics; rhandle=@riccati; % set up handles for riccati eqn and dynamic eqn
thresh=1.e-3; % termination threshold
sty={'' '--' '-.' ':'}; % set up line styles for plotting
%for computing values of control signal(s)
dT = 0.01;
To = 0; Tf = 10;
tspan = To:dT:Tf;
old_u=zeros(1,length(tspan));
xo=[0 0 0 0]'; % set up IC - the first element will change depending on which level of pathogen attack is used
for ii = 1:length(IC) % loop round the different pathogen attack scenarios
  term=100; % re-set termination initial value for each run
  solx=[];solp=[]; % erase previous values of x(t) and P(t)
  xo(1)=IC(ii);solx=xo; % First iteration uses fixed value of xo.
  cnt=0; % counter
  while term>thresh
      solp=ode23s(rhandle,[0 Tf],po,options,sys_name,solx,Tf); % solve riccati ode
      solx=ode23s(dhandle,[0 Tf],xo,options,sys_name,solx,solp,Tf); % solve dynamics ode
      [u,J,Jf]=evalcontrol(tspan,sys_name,solx,solp,Tf); % compute control signal
      term=(norm(u-old_u))/Tf; % compute termination test value
      old_u=u; % over-write old controls with new
      disp(term); cnt=cnt+1;
  end
  disp(['No. Iterations = ' int2str(cnt)])
subplot(3,2,1);plot(solx.x,solx.y(1,:),sty{ii},'linewidth',1);axis([0 10 0 5]);ylabel('pathogen conc.');grid on;hold on
subplot(3,2,2);plot(solx.x,solx.y(2,:)+sys.xe(2),sty{ii},'linewidth',1);axis([0 10 0 10]);ylabel('plasma conc.');grid on;hold on
subplot(3,2,3);plot(solx.x,solx.y(3,:)+sys.xe(3),sty{ii},'linewidth',1);axis([0 10 0 5]);ylabel('antibody conc.');grid on;hold on
subplot(3,2,4);plot(solx.x,solx.y(4,:),sty{ii},'linewidth',1);axis([0 10 0 1]);ylabel('organ health');grid on;hold on
subplot(3,2,5);plot(tspan,u,sty{ii},'linewidth',1);ylabel('controls');grid on;hold on
subplot(3,2,6);plot(tspan,J,sty{ii},'linewidth',1);ylabel('cost');grid on;hold on
end
legend('sub-clinical','clinical','chronic','lethal','location','best');

I also have f_of_x.m, which I will include below.

function dxdt=f_of_x(t,x,dT,u)
a11=1;a12=1;a23=1;a31=1;a42=1;b2=1;b3=1;
b1=-1;b4=-1;
a22=3;
a32=1.5;
a33=0.5;a41=0.5;
x2bar=2;
n=floor(t/dT)+1;
dxdt=zeros(4,1);
dxdt(1)=(a11-a12*x(3))*x(1)+b1*u(n,1);
dxdt(2)=a21(x(4))*a22*x(1)*x(3)-a23*(x(2)-x2bar)+b2*u(n,2);
dxdt(3)=a31*x(2)-(a32+a33*x(1))*x(3)+b3*u(n,3);
dxdt(4)=a41*x(1)-a42*x(4)+b4*u(n,4);
function y=a21(x)
if x>=0 & x<0.5
  y=cos(pi*x);
elseif x>=0.5
  y=0;
else
  error('negative values not allowed')
end

[Walter Roberson]

Looks like we still need imresp2

[chromaclouds]

oops sorry, imresp2 is imresp. I forgot to change the name of the file. everything that is 'imresp2' is actually 'imresp' (in evalcontrol and ACS420_3).

Answer 1

[Walter Roberson]

In imresp you have

sys.b = [b1 0 0 0; 0 b2 0 0]';

that forces sys.b to be 4 x 2. When you work with in in evalcontrol in

      u(:,ii)=-sys.r\sys.b'*reshape(deval(solp,Tf-t),length(x),length(x))*x;

then that forces there to be two rows of output. The various 1 x 4 and 4 x 4 cancel out effectively to make the number of columns of sys.b the control over the number of rows of result.

So the right hand side ends up being 1 x 2, being stored into column #ii of u. And that is a problem because you initialized

u=zeros(1,length(tspan));

so u only has one row, but you are trying to store two rows.

If you just change to

u=zeros(2,length(tspan));

then the code will proceed on to some plots.

[chromaclouds]

that did it! thank you very much!