Matlab fatal error after/while executing "grampc_run_Cmex"

[art ber]

Hello,

I tried to run the following code:
startMPC.m

function [vec,grampc] = startMPC()
  %
  % This file is part of GRAMPC.
  %
  % GRAMPC - a gradient-based MPC software for real-time applications
  %
  % Copyright (C) 2014 by Bartosz Kaepernick, Knut Graichen, Tilman Utz
  % Developed at the Institute of Measurement, Control, and
  % Microtechnology, University of Ulm. All rights reserved.
  %
  % GRAMPC is free software: you can redistribute it and/or modify
  % it under the terms of the GNU Lesser General Public License as 
  % published by the Free Software Foundation, either version 3 of 
  % the License, or (at your option) any later version.
  %
  % GRAMPC is distributed in the hope that it will be useful,
  % but WITHOUT ANY WARRANTY; without even the implied warranty of
  % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  % GNU Lesser General Public License for more details.
  %
  % You should have received a copy of the GNU Lesser General Public 
  % License along with GRAMPC. If not, see <http://www.gnu.org/licenses/>.
  %
  % File: startMPC.m
  % Authors: Bartosz Kaepernick, Knut Graichen, Tilman Utz
  % Date: February 2014
  % Version: v1.0
  %
  % MATLAB/Simulink file for running GRAMPC example.
  %

close all; clear all; clc;
  %% Compilation
  PROBFCT = 'probfct_test.c';
  COMPILE = 1;
  DEBUG   = 0;
  VERBOSE = 0;
  CLEAN   = 0;

  if COMPILE || CLEAN
    if exist('make','file') ~= 2
      curPath = cd;
      cd('..');
      addpath(cd);
      cd(curPath);
    end
  end

  if CLEAN
    make clean;
    vec    = [];
    grampc = [];
    return;
  end

  if COMPILE
    if VERBOSE && DEBUG
      make(PROBFCT,'verbose','debug');
    elseif VERBOSE
      make(PROBFCT,'verbose');
    elseif DEBUG
      make(PROBFCT,'debug');
    else
      make(PROBFCT);
    end
  end

  %% Initialization
  grampc = grampc_init_Cmex();

  %% user options (all options are optional)
  % grampc_setopt_Cmex(grampc,'MaxIter',5);
  % grampc_setopt_Cmex(grampc,'ShiftControl','on');
  % grampc_setopt_Cmex(grampc,'ScaleProblem','on');
  grampc_setopt_Cmex(grampc,'CostIntegrator','simpson');
  grampc_setopt_Cmex(grampc,'Integrator','euler');
  % grampc_setopt_Cmex(grampc,'IntegratorRelTol',1e-3);
  % grampc_setopt_Cmex(grampc,'IntegratorAbsTol',1e-3);
  grampc_setopt_Cmex(grampc,'LineSearchType','adaptive');
  % grampc_setopt_Cmex(grampc,'LineSearchMax',2.0);
  % grampc_setopt_Cmex(grampc,'LineSearchMin',1e-3);
  % grampc_setopt_Cmex(grampc,'LineSearchInit',0.1);
  % grampc_setopt_Cmex(grampc,'LineSearchIntervalFactor',0.5);
  % grampc_setopt_Cmex(grampc,'LineSearchAdaptFactor',2.0);
  % grampc_setopt_Cmex(grampc,'LineSearchIntervalTol',0.2);
  grampc_setopt_Cmex(grampc,'JacobianX','sysjacxadj');
  grampc_setopt_Cmex(grampc,'JacobianU','sysjacuadj');
  grampc_setopt_Cmex(grampc,'IntegralCost','on');
  grampc_setopt_Cmex(grampc,'FinalCost','on');

  %% User parameter
  % mandatory paramter 
  user.param.x0     = [0.0, 0.0];
  user.param.u0     = [0.0, 0.0];
  user.param.xdes   = [-66.0 231.1];
  user.param.udes   = [-67.6909, 36.5465];
  user.param.Thor   = 1e-5;
  user.param.dt     = 1e-5;
  % optional parameter  
  user.param.umax   =  2/3 * 250 * [1.0, 1.0];
  user.param.umin   = -2/3 * 250 * [1.0, 1.0];
  user.param.Nhor   = 1;
  user.param.NpCost = 6;
  user.param.pCost  = [1.0, 1.0,...
                           1.0, 1.0,...
                       0.0, 0.0];

  grampc_setparam_Cmex(grampc,'xk',user.param.x0);
  grampc_setparam_Cmex(grampc,'u0',user.param.u0);
  grampc_setparam_Cmex(grampc,'xdes',user.param.xdes);
  grampc_setparam_Cmex(grampc,'udes',user.param.udes);
  grampc_setparam_Cmex(grampc,'umax',user.param.umax);
  grampc_setparam_Cmex(grampc,'umin',user.param.umin);
  grampc_setparam_Cmex(grampc,'Thor',user.param.Thor);
  grampc_setparam_Cmex(grampc,'Nhor',user.param.Nhor);
  grampc_setparam_Cmex(grampc,'dt',user.param.dt);
  grampc_setparam_Cmex(grampc,'NpCost',user.param.NpCost);
  grampc_setparam_Cmex(grampc,'pCost',user.param.pCost);

  %% MPC loop

  Tsim = 7e-4;

  kk = 1;
  vec.t = 0:grampc.param.dt:Tsim;  
  vec.u = zeros(grampc.param.Nu,length(vec.t));
  vec.x = zeros(grampc.param.Nx,length(vec.t));
  vec.x(:,1) = grampc.param.xk;
  vec.J = zeros(1,length(vec.t));
  vec.CPUtime = zeros(1,length(vec.t)-1);

  while (vec.t(kk)<Tsim)

    tic;
    [xnext,unext,J] = grampc_run_Cmex(grampc);
    grampc_setparam_Cmex(grampc,'xk',xnext);
    vec.CPUtime(kk) = toc;

    vec.x(:,kk+1) = xnext;
    vec.u(:,kk+1) = unext;
    vec.J(kk+1)   = J;

    kk = kk + 1;

  end

  figure(1);

  subplot(2,1,1);
  plot(vec.t,vec.x);
  grid on;

  subplot(2,1,2);
  plot(vec.t, vec.u);
  grid on;

% ******* END OF FUNCTION *******
end

probfct_test.c

/*
 *
 * This file is part of GRAMPC.
 *
 * GRAMPC - a gradient-based MPC software for real-time applications
 *
 * Copyright (C) 2014 by Bartosz Kaepernick, Knut Graichen, Tilman Utz
 * Developed at the Institute of Measurement, Control, and
 * Microtechnology, University of Ulm. All rights reserved.
 *
 * GRAMPC is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as 
 * published by the Free Software Foundation, either version 3 of 
 * the License, or (at your option) any later version.
 *
 * GRAMPC is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public 
 * License along with GRAMPC. If not, see <http://www.gnu.org/licenses/>.
 *
 */

/*
 * 
 * File: probfct.c
 * Authors: Bartosz Kaepernick, Knut Graichen, Tilman Utz
 * Date: February 2014
 * Version: v1.0
 *
 * Problem function for grampc toolbox.
 *
 *
 * This probfct file describes the dynamics, the cost function and the corresponding
 * derivatives of an overhead crane with 6 states and 2 controls. For a more
 * detailed model, for instance, see
 *
 * B. Käpernick and K. Graichen, "Model Predictive Control of an Overhead Crane
 * Using Constraint Substitution", Proceedings of the American Control Conference,
 * Washington D.C. (USA), 2013.
 *
 *
 * This probfct file provides an interface to GRAMPC. The underlying
 * optimal control problem (OCP) of the model predictive control (MPC) formulation
 * has the following structure
 *                                  _T
 *                                 /
 *      min    J(u,xk) = V(T,x(T)) + / L(t,x(t),u(t)) dt
 *      u(.)                    _/
 *                             0
 *             .
 *      s.t.   x(t) = f(tk+t,x(t),u(t)), x(0) = xk
 *
 *             Ul <= u(t) <= Uu, t in [0,T]
 *
 * with states x(t), constrained controls u(t) and the fixed prediction horizon T.
 * The functions V(t,x), L(t,x,u) and f(t,x,u) denote the terminal and integral
 * cost and the systems dynamics. Note that no terminal conditions for the states
 * are included in the problem formulation.
 *
 * The necessary optimality conditions can then be derived by means of the
 * Hamiltonian
 *
 * The function interfaces below have the following meaning (adj denotes the costates):
 *
 * sysfct:     f(t,x,u)
 *
 * sysjacxadj: df(t,x,u)/dx' * adj
 *
 * sysjacuadj: df(t,x,u)/du' * adj
 *
 * sysjacx:    df(t,x,u)/dx
 *
 * sysjacu:    df(t,x,u)/du
 *
 * icostfct:   L(t,x,u)
 *
 * icostjacx:  dL(t,x,u)/dx
 *
 * icostjacu:  dL(t,x,u)/du
 *
 * fcostfct:   V(t,x)
 *
 * fcostjacx:  dV(t,x)/dx
 *
 */

#include "probfct.h"

#define a11 -111.815132314573
#define a12  1059.45323629074
#define a21 -662.453528997459
#define a22 -88.4173297966401
#define b1   3727.17107715244
#define b2   2947.24432655467
#define w   -130861.114561555

void sysdim(typeInt *Nx, typeInt *Nu)
{
  *Nx = 2;   // Number of states
  *Nu = 2;   // Number of control inputs
}

void sysfct(typeRNum *out, typeRNum t, typeRNum *x, typeRNum *u, typeRNum *pSys)
{
    out[0] = a11*x[0] + a12*x[1] + b1*u[0];      // f0 := x0'= a11*x0 + a12*x1 + b1*u0
    out[1] = a21*x[0] + a22*x[1] + b2*u[1] + w;  // f1 := x1'= a21*x0 + a22*x1 + b2*u1 + w
}

void sysjacxadj(typeRNum *out, typeRNum t, typeRNum *x, typeRNum *adj, typeRNum *u, typeRNum *pSys)
{
    out[0] = a11*adj[0] + a21*adj[1];    // (df0/dx0)*adj0 + (df1/dx0)*adj1
    out[1] = a12*adj[0] + a22*adj[1];    // (df0/dx1)*adj0 + (df1/dx1)*adj1
}

void sysjacuadj(typeRNum *out, typeRNum t, typeRNum *x, typeRNum *adj, typeRNum *u, typeRNum *pSys)
{
    out[0] = b1*adj[0];    // (df0/du0)*adj0 + (df1/du0)*adj1
    out[1] = b2*adj[1];    // (df0/du1)*adj0 + (df1/du1)*adj1
}

/*
 * Matrix must be entered rowwise, i.e. A=[a11,a12;a21,a22] becomes
 * out[0]=a11; out[1]=a12; out[2]=a21; out[3]=a22;
 * Note that Matlab/Fortran does it columnwise !
 */
void sysjacx(typeRNum *out, typeRNum t, typeRNum *x, typeRNum *u, typeRNum *pSys)
{
    out[0] = a11;   // df0/dx0
    out[1] = a21;   // df1/dx0
    out[2] = a12;   // df0/dx1
    out[3] = a22;   // df1/dx1
}

/*
 * Matrix must be entered rowwise, i.e. A=[a11,a12;a21,a22] becomes
 * out[0]=a11; out[1]=a12; out[2]=a21; out[3]=a22;
 * Note that Matlab/Fortran does it columnwise !
 */
void sysjacu(typeRNum *out, typeRNum t, typeRNum *x, typeRNum *u, typeRNum *pSys)
{
    out[0] = b1;       // df0/du0
    out[1] = 0.0;      // df1/du0
    out[2] = 0.0;      // df0/du1
    out[3] = b2;       // df1/du1
}

void fcostfct(typeRNum *out, typeRNum t,  typeRNum *x, typeRNum *xdes, typeRNum *pCost)
{                                                              // V = (Delta x )^T * P * (Delta x)
    out[0] = ( x[0] - xdes[0] )*( x[0] - xdes[0] )*pCost[0]    //   = (Delta x0)^2 * P0
           + ( x[1] - xdes[1] )*( x[1] - xdes[1] )*pCost[1];   //   + (Delta x1)^2 * P1
}

void fcostjacx(typeRNum *out, typeRNum t, typeRNum *x, typeRNum *xdes, typeRNum *pCost)
{
    out[0] = 2*( x[0] - xdes[0] )*pCost[0];   // dV/d(Delta x0)
    out[1] = 2*( x[1] - xdes[1] )*pCost[1];   // dV/d(Delta x1)
}

void icostfct(typeRNum *out, typeRNum t, typeRNum *x, typeRNum *u, typeRNum *xdes, typeRNum *udes, typeRNum *pCost)
{                                                              // L = (Delta x )^T * Q * (Delta x) + (Delta u)^T * R * (Delta u)
    out[0] = ( x[0] - xdes[0] )*( x[0] - xdes[0] )*pCost[2]    //   = (Delta x0)^2 * Q0
           + ( x[1] - xdes[1] )*( x[1] - xdes[1] )*pCost[3]    //   + (Delta x1)^2 * Q1
           + ( u[0] - udes[0] )*( u[0] - udes[0] )*pCost[4]    //   + (Delta x1)^2 * R0
           + ( u[1] - udes[1] )*( u[1] - udes[1] )*pCost[5];   //   + (Delta x1)^2 * R1
}

void icostjacx(typeRNum *out, typeRNum t, typeRNum *x, typeRNum *u, typeRNum *xdes, typeRNum *udes, typeRNum *pCost)
{
    out[0] = 2*( x[0] - xdes[0] )*pCost[2];    // dL/d(Delta x0)
    out[1] = 2*( x[1] - xdes[1] )*pCost[3];    // dL/d(Delta x1)
}

void icostjacu(typeRNum *out, typeRNum t, typeRNum *x, typeRNum *u, typeRNum *xdes, typeRNum *udes, typeRNum *pCost)
{
    out[0] = 2*( u[0] - udes[0] )*pCost[4];     // dL/d(Delta u0)
    out[1] = 2*( u[1] - udes[1] )*pCost[5];     // dL/d(Delta u1)
}

After compiling when it comes to the simulation I get a fatal error.

Can anyone help me?

[Tobias Englert]

Hello,

Nhor must be greater than 1.

Best regards

[art ber]

Oh no! Thank you very much!