Why fmincon did not work when it is provided with gradient

Hi all; I have a problem with two variables. If the objective function is quadratic and subjected to boundary constraints and dynamic constraints which are consider as a linear constraints. I know the objective function is a function that returns a scalar my question is if i provide the gradient of the objective function which should be a vector why i get this message error.

Warning: Input arguments must be scalar. 
> In objfunTetsuro at 27
  In @(kc)objfunTetsuro(kc)
  In fmincon at 599
  In main at 26
???  In an assignment  A(:) = B, the number of elements in A and B
must be the same.

Error in ==> fmincon at 599
      [initVals.f,initVals.g(:)] = feval(funfcn{3},X,varargin{:});

Error in ==> main at 26
[kc_opt,fopt,flag,output]=fmincon(@(kc)objfunTetsuro(kc),x0,[],[],[],[],lb,ub,[],...

Caused by:
    Failure in initial user-supplied objective function evaluation.
    FMINCON cannot continue.

where the objective function is

function [val,GObj] = objfunTetsuro(k_c)

load IRI_737b
zdot = calc_zdot(road_x, road_z, 10); 
zt=zdot.u;

x00=[0,0,0,0];

opt = odeset('RelTol', 1e-2, 'AbsTol', 1e-3);
[t02,y] = ode23(@(t,x)f(t,x,0, k_c, @(t)lookup_u(zdot,t)), [0 2], x00,opt);
x_frame=y';

val = 0;
R = Weight();
x_frame_zs = zeros(5,length(t02));

xx=size(t02,1);

for i=1:(xx-1)
    ddx = ff(t02(i),x_frame(1:4,i), 0, k_c, zt(i));
    x_frame_zs(:,i) = [x_frame(:,i); ddx(4)];
    val = val+ (t02(i+1)-t02(i))*x_frame_zs(:,i)'*R*x_frame_zs(:,i);
end

   if nargout > 1
      R_gradient = Weight_gradient(k_c);
      GObj = zeros(t02,1);  % line 27 which is the line error

     for i=1:(t02-1)
         GObj((i-1)*4+1:i*4) = 2*(t(i+1)-t(i))*R_gradient*x_frame(:,i);
     end 
   end
end

and fmincon is

[kc_opt,fopt,flag,output]=fmincon(@(kc)objfunTetsuro(kc),x0,[],[],[],[],lb,ub,[],...  % line error 28
    optimset('GradObj','on','Display','iter','Algorithm','Interior-Point',...
             'DiffMinChange',1e-3,'DiffMaxChange',1));

Answer 1

[Walter Roberson]

You have

[t02,y] = ode23(@(t,x)f(t,x,0, k_c, @(t)lookup_u(zdot,t)), [0 2], x00,opt);

That is going to result in t02 being a column vector of times within your timespan [0 2].

Then later you have

      GObj = zeros(t02,1);  % line 27 which is the line error

as if t02 is a scalar integer that is a size. But it is not a scalar and it is not an integer.

Furthermore right after that you have

      for i=1:(t02-1)

where again the context would require t02 to be an integer scalar. But it isn't.

Looking at your code, I would suggest it should be

     GObj = zeros(xx,1);  % line 27 which is the line error

     for i=1:(xx-1)

[Muna Shehan]

Thanks alot for your reply; I edit t02 but i still get a massage error. Kindly can you explain why i get this error

???  In an assignment  A(:) = B, the number of elements in A and B
must be the same.

Error in ==> fmincon at 599
      [initVals.f,initVals.g(:)] = feval(funfcn{3},X,varargin{:});

Error in ==> main at 26
[kc_opt,fopt,flag,output]=fmincon(@(kc)objfunTetsuro(kc),x0,[],[],[],[],lb,ub,[],...

Caused by:
    Failure in initial user-supplied objective function evaluation.
    FMINCON cannot continue.

[Walter Roberson]

I would need your complete source and the content of IRI_737b to test.

[Muna Shehan]

Thank you Walter for your reply. This is fmincon with lb, ub and x0

x0=[1000,550];
lb=[1000,500];
ub=[20000,4000];
[kc_opt,fopt,flag,output]=fmincon(@(kc)objfunTetsuro(kc),x0,[],[],[],[],lb,ub,[],...
    optimset('GradObj','on','Display','iter','Algorithm','Interior-Point','DiffMinChange',1e-3,'DiffMaxChange',1));

where objfunTetsuro returns the objective function with gradient ( objfunTetsuro is saved as objfunTetsuro.m in the path)

function [val,GObj] = objfunTetsuro(k_c)

load IRI_737b
zdot = calc_zdot(road_x, road_z, 10); 
zt=zdot.u;

x00=[0,0,0,0];  

opt = odeset('RelTol', 1e-2, 'AbsTol', 1e-3);
[t02,y] = ode23(@(t,x)f(t,x,0, k_c, @(t)lookup_u(zdot,t)), [0 2], x00,opt);
x_frame=y';

val = 0;
R = Weight();
x_frame_zs = zeros(5,length(t02));

xx=size(t02,1);

for i=1:(xx-1)%(length(t02)-1)
    ddx = ff(t02(i),x_frame(1:4,i), 0, k_c, zt(i));
    x_frame_zs(:,i) = [x_frame(:,i); ddx(4)];
    val = val+ (t02(i+1)-t02(i))*x_frame_zs(:,i)'*R*x_frame_zs(:,i);
end

   if nargout > 1
      R_gradient = Weight_gradient(k_c);
      GObj = zeros(xx,1);
   for i=1:(xx-1)
      GObj((i-1)*4+1:i*4) = 2*(t02(i+1)-t02(i))*R_gradient*x_frame(:,i);
  end 
  end
  end

f(t,x,0, k_c, @(t)lookup_u(zdot,t) is a handle function used to solve the ODE for the dynamic system which are consider as a linear constraints because my problem is similar to what Mathwork support describe in the link below and I try to implement the first approache. http://www.mathworks.com/matlabcentral/answers/101883-how-do-i-estimate-or-optimize-the-parameters-of-my-ode-system-in-matlab-8-1-r2013a

function dx = f(t, x, ~, kc,z)
param.ms=240;
param.mus=36;
param.kus=16e4;
param.ct = 0 ;          

u=0;
J = jacobian(t, kc, param);
A = J(1:4,1:4);
B = J(1:4, 5);
b_ramp = [-1; param.ct/param.mus; 0; 0];

if ~isempty(z)   
    dx= A* x + B * u + b_ramp * z(t); 
else
    dx = A*x + B*u ;
end
end

and

function u = lookup_u(tu, t)
i = find(tu.t <= t, 1, 'last');
[m,n] = size(tu.u);
if (m==1) || (n==1)
    if i<length(tu.t)
        u = (tu.u(i)*(tu.t(i+1) -t) + tu.u(i+1)*(t-tu.t(i)))/(tu.t(i+1)-tu.t(i));
    else
        u = tu.u(i); 
    end
else
    if i < length(tu.t)
        u = (tu.u(:,i)*(tu.t(i+1) -t) + ...
            tu.u(:,i+1)*(t-tu.t(i)))/(tu.t(i+1)-tu.t(i));
    else
        u = tu.u(:,i);
    end
end
end

and jacobian is

function J = jacobian(~, kc, param)
ks=kc(1);
cs=kc(2);
A = [ 0 1 0 0;
    -param.kus/param.mus -cs/param.mus ks/param.mus cs/param.mus;
    0 -1 0 1;
    0 cs/param.ms -ks/param.ms -cs/param.ms];
B = [0 -1/param.mus 0 1/param.ms]'; 
J = [A, B];
end

and ff(t, x, ~, kc,z) function is used to estimate the states variable for each time step (z is a the vertical velocity that disturb the system).

function ddx = ff(t, x, ~, kc,z)
param.ms=240;
param.mus=36;
param.kus=16e4;
param.ct = 0 ;          
u=0;
J = jacobian(t, kc, param);
A = J(1:4,1:4);
B = J(1:4, 5);
b_ramp = [-1; param.ct/param.mus; 0; 0];
if ~isempty(z)
  ddx= A* x + B * u + b_ramp * z; 
else
  ddx = A*x + B*u ;
end
end

and Weight() is

function R = Weight()
r1 = 1e+5;
r2 = 0.005;
R = diag([r1, 0, 0, 0, r2]);
end

and Weight_gradient ()

function R_gradient = Weight_gradient(xd)
param.ms=240;
param.mus=36;
param.kus=16e4;
param.ct = 0 ;          

r1 = 1e+5;
r2 = 0.005;
ks = xd(1);
cs = xd(2);
R_gradient = [
    r1 0 0 0 ;...
    0 r2*(cs/param.ms)^2 -r2*(cs/param.ms)*(ks/param.ms) -r2*(cs/param.ms)^2 ;...
    0 -r2*(cs/param.ms)*(ks/param.ms) r2*(ks/param.ms)^2 r2*(cs/param.ms)*(ks/param.ms) ;...
    0 -r2*(cs/param.ms)^2 r2*(cs/param.ms)*(ks/param.ms) r2*(cs/param.ms)^2 ];
end

and calc_zdot

function out = calc_zdot(x, z, v)
% Calculate dz/dt from z-t data.
n = length(x);
t = zeros(length(x),1);
zdot = zeros(length(x)-1, 1);
for i=2:n
    dist = sqrt((z(i) - z(i-1))^2 + (x(i)-x(i-1))^2);
    t(i) = t(i-1) + dist/v;
    zdot(i-1) = (z(i) - z(i-1))/(t(i) - t(i-1));
end
out = struct('t', t, 'u', zdot);
end

These are the all functions that I used in my code and I attached the road excitation file (IRI_737d). and its the only way that I know to attach the source code, if there is another way please advise me. Thanks again for your time.