Choise of prediction horizon

[hcl734]

Hi all,

I have a question regarding the choice of prediction horizon.
In my understanding a bigger prediction horizon should stabilize the control of the system in MPC.
However when I choose the prediction horizon (N Ts) for my example below to high the results of a setpoint tracking simulation deteriorate with growing prediction horizon.
I attached the results of the simulation for prediction horizons of 60, 90 and 100 seconds to demonstrate my problem. I figured out then that if I raise the number of iterations according to the length of the predicition horizon (N Ts) the results get equivalent for all horizon length

user.opt.Nhor        = N*Ts;          % Number of steps for the system integration
user.opt.MaxGradIter = N*Ts;           % Maximum number of gradient iterations
user.opt.MaxMultIter = N*Ts;           % Maximum number of augmented Lagrangian iterations
user.opt.ConvergenceCheck = 'on';

When I choose a prediction horizon higher than exactly 112 my simulation terminates becaue of NaN values although my problemfct.c does not define divisions by zero even for unphysical values of states and controls.
Does somebody know how those NaN values could be explained?

Also if I use grampc_estim_penmin_Cmex for a horizon bigger than 112 it gives the error:
PenaltyMin: Invalid value for option.
Error using grampc_estim_penmin_Cmex
Invalid setting of option or parameter

I am clearly missing something here. Any hint is appreciated.
Best regards

startMPC.m

function [vec,grampc,figNr] = startMPC(figNr,compile,varargin)

%% Check input arguments
if nargin < 1 || isempty(figNr)
    figNr = 1;
end
if nargin < 2 || isempty(compile)
    compile = 0;
end

%% Parameters
% path to grampc root
grampc_root_path = '../GRAMPC_v2.1/';
addpath([grampc_root_path 'matlab/mfiles']);
% name of problem function
probfct = 'probfct_tank.c';

% plot predicted trajectories
PLOT_PRED = 1;
% plot solution trajectories
PLOT_TRAJ = 1;
% plot optimization statistics
PLOT_STAT = 1;
% update plots after N steps
PLOT_STEPS = 10;
% pause after each plot
PLOT_PAUSE = 0;

% problem-specific plot
PLOT = 0;

% Options for the reference simulation
odeopt =  [];%odeset('RelTol',1e-6,'AbsTol',1e-8);

%% Compilation
% compile toolbox
if compile > 1 || ~exist([grampc_root_path 'matlab/bin'], 'dir')
    grampc_make_toolbox(grampc_root_path, varargin{:});
end
% compile problem
if compile > 0 || ~exist('+CmexFiles', 'dir')
    grampc_make_probfct(grampc_root_path, probfct, varargin{:});
end

%% Initialization
% init GRAMPC and print options and parameters
[grampc,Tsim] = initData;
CmexFiles.grampc_printopt_Cmex(grampc);
CmexFiles.grampc_printparam_Cmex(grampc);

% init solution structure
vec = grampc_init_struct_sol(grampc, Tsim);

% init plots and store figure handles
if PLOT_PRED
    phpP = grampc_init_plot_pred(grampc,figNr);
    figNr = figNr+1;
end
if PLOT_TRAJ
    phpT = grampc_init_plot_sim(vec,figNr);
    figNr = figNr+1;
end
if PLOT_STAT
    phpS = grampc_init_plot_stat(vec,grampc,figNr);
    figNr = figNr+1;
end

%% MPC loop
i = 1;
while 1
    % set current time and current state
    grampc = CmexFiles.grampc_setparam_Cmex(grampc,'t0',vec.t(i));
    grampc = CmexFiles.grampc_setparam_Cmex(grampc,'x0',vec.x(:,i));

%     % Führungsgrößensprung
    if i > 200
        grampc = CmexFiles.grampc_setparam_Cmex(grampc,'xdes', [150, 130, 125, 5000, 5000, 5000]);
    end
    % run MPC and save results
    [grampc,vec.CPUtime(i)] = CmexFiles.grampc_run_Cmex(grampc);
    vec = grampc_update_struct_sol(grampc, vec, i);

    % print solver status
    printed = CmexFiles.grampc_printstatus_Cmex(grampc.sol.status,'Error');
    if printed
        fprintf('at simulation time %f.\n --------\n', vec.t(i));    
    end

    % check for end of simulation
    if i+1 > length(vec.t)
        break;
    end

    % simulate system
    [~,xtemp] = ode45(@CmexFiles.grampc_ffct_Cmex,vec.t(i)+[0 double(grampc.param.dt)],vec.x(:,i),odeopt,vec.u(:,i),vec.p(:,i),grampc.userparam);
    vec.x(:,i+1) = xtemp(end,:);

    % evaluate time-dependent constraints 
    % to obtain h(x,u,p) instead of max(0,h(x,u,p))
    vec.constr(:,i) = CmexFiles.grampc_ghfct_Cmex(vec.t(i), vec.x(:,i), vec.u(:,i), vec.p(:,i), grampc.userparam);

    % update iteration counter
    i = i + 1;

    % plot data
    if mod(i,PLOT_STEPS) == 0 || i == length(vec.t)
        if PLOT_PRED
            grampc_update_plot_pred(grampc,phpP);
        end
        if PLOT_TRAJ
            grampc_update_plot_sim(vec,phpT);
        end
        if PLOT_STAT
            grampc_update_plot_stat(vec,grampc,phpS);
        end
        drawnow
        if PLOT_PAUSE
            pause;
        end
    end
end
end

initData.m

function [grampc,Tsim] = initData()

%% Allgemeine Parameter
Ts = 0.5; % Sampling Zeit
N = 224; % control intervals
Tf = 600; % Simulationsdauer

%% Parameter definition
% Initial values and setpoints of the states
user.param.x0    = [111, 100, 95, 3.0, 3.4, 6];
user.param.xdes  = [111, 100, 95, 5000, 5000, 5000]; %letzte drei Einträge sind dummies, da Dimension von xdes = Nx sein muss

% Initial values, setpoints and limits of the inputs
user.param.u0    = [0.0, 0.0, 0.0];
user.param.udes  = [0.0, 0.0, 0.0];
user.param.umax  = [9.9, 9.9, 2]; % ersetzen Beschränkungen für Stellgrößen der MPC (delta_)
user.param.umin  = [-9.9, -9.9, -2];

% Time variables
user.param.Thor  = N*Ts;           % Prediction horizon
user.param.dt    = Ts;             % Sampling time
user.param.t0    = 0.0;             % time at the current sampling step

%% Option definition
% Basic algorithmic options
user.opt.Nhor        = N*Ts;          % Number of steps for the system integration
user.opt.MaxGradIter = N*Ts;           % Maximum number of gradient iterations
user.opt.MaxMultIter = N*Ts;           % Maximum number of augmented Lagrangian iterations
user.opt.ConvergenceCheck = 'on';

% Constraint tolerances
user.opt.ConstraintsAbsTol = 1e-3*[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
% user.opt.ConstraintsAbsTol = [0.1 104.5];

% optional settings for a better performance
% user.opt.PenaltyMin = 10000; % Comment line 83 (grampc = CmexFiles.grampc_estim_penmin_Cmex(grampc,1);) to use this option value
% user.opt.AugLagUpdateGradientRelTol = 1e0;
% user.opt.LineSearchMax = 1e1;

% System integration
user.opt.TerminalCost = 'on';

%% Userparameter definition
% e.g. system parameters or weights for the cost function
p(1) = 5.977e-01;
p(2) = 7.286e-01;
p(3) = 5.2344e-04;
p(4) = 4.7e-03;
p(5) = 5.1991e-04;
p(6) = 6.374e-01;
p(7) = 9.81e+03;

Wh1 = 1;
Wh2 = 1;
Wh3 = 1;
WhN1 = 1;
WhN2 = 1;
WhN3 = 1;
Wdelta = 10;
Wdeltap1 = Wdelta;
Wdeltap2 = Wdelta;
Wdeltav12 = Wdelta;

pCost = [Wh1, Wh2, Wh3, WhN1, WhN2, WhN3, Wdeltap1, Wdeltap2, Wdeltav12]; %size(pCost) = 9

pConstr_x_up = [130, 130, 130, 4.95, 4.95, 4.95]; %size(pConstr_x_up) = 6
pOffset_x_up = [160, 160, 160, 4.95, 4.95, 4.95];
pScale_x_up = pConstr_x_up;
pConstr_x_down = -pConstr_x_up;%size(pConstr_x_up) = 6
pOffset_x_down = pOffset_x_up;
pScale_x_down = pConstr_x_up;
pConstr_x = [pConstr_x_up pConstr_x_down];

pConstr_u_up = [9.9, 9.9, 2];
pOffset_u_up = [5000, 5000, 5000]; % nicht implementiert!
pScale_u_up = pConstr_u_up;
pConstr_u_down = [-9.9, -9.9, -2];
pOffset_u_down = [5000, 5000, 5000]; % nicht implementiert!
pScale_u_down = pConstr_u_up;
pConstr_u = [pConstr_u_up pConstr_u_down];

pOffset = [pOffset_x_up, pOffset_x_down, pOffset_u_up, pOffset_u_down]; % [18x1]
pScale = [pScale_x_up, pScale_x_down, pScale_u_up, pScale_u_down]; % [18x1]

userparam = [p pCost pConstr_x pConstr_u pOffset pScale];

%% Grampc initialization
grampc = CmexFiles.grampc_init_Cmex(userparam);

%% Update grampc struct while ensuring correct data types
grampc = grampc_update_struct_grampc(grampc,user);

%% Estimate and set PenaltyMin
grampc = CmexFiles.grampc_estim_penmin_Cmex(grampc,1);
grampc.opt.PenaltyMin
%% Simulation time
Tsim = Tf;

probfct.c

#include "probfct.h"

#define NPSYS 7 // size(p), s.u.
#define NPCOST 9 // size(p), s.u.
#define NPCONSTR 18
#define NPOFFSET 18

 /* square macro */
#define POW2(a) ((a)*(a))
#define SIGNAP(x)                   tanh(10000*x)
#define ABSAP(x)                sqrt(pow(x,2)+1e-10)
#define TANH(x)                     tanh(x)
#define SQRT(x)           ((x) > 0 ? sqrt(x): 1e-16)
#define ALMOSTZERO(x)           ((x) != 0 ? (x): 1e-16)

/** OCP dimensions: states (Nx), controls (Nu), parameters (Np), equalities (Ng),
    inequalities (Nh), terminal equalities (NgT), terminal inequalities (NhT) **/
void ocp_dim(typeInt *Nx, typeInt *Nu, typeInt *Np, typeInt *Ng, typeInt *Nh, typeInt *NgT, typeInt *NhT, typeUSERPARAM *userparam)
{
    *Nx = 6; // h1=x[0] h2=x[1] h3=x[2] p1=x[3] p2=x[4] v12=x[5]
    *Nu = 3; // delta_p1=u[0] delta_p2=u[1] delta_v12=u[2]
    *Np = 0; //
    *Nh = 18; //
    *Ng = 0;
    *NgT = 0;
    *NhT = 0;
}

/** System function f(t,x,u,p,userparam)
    ------------------------------------ **/
void ffct(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
    // typeRNum ist double
    // ctypeRNum ist constant double
    ctypeRNum *pSys = (typeRNum*)userparam;
    ctypeRNum f3 = 2;
    // C arrays werden ab 0 indiziert --> p(1) = pSys[0]...; p(7) = pSys[6]; p(6) = pSys[5]
    out[0] = pSys[0]*x[3]-x[5]*pSys[2]*SIGNAP(x[0]-x[1])*SQRT(f3*pSys[6]*ABSAP(x[0]-x[1]));
    out[1] = x[5]*pSys[2]*SIGNAP(x[0]-x[1])*SQRT(f3*pSys[6]*ABSAP(x[0]-x[1]))-pSys[3]*SIGNAP(x[1]-x[2])*SQRT(f3*pSys[6]*ABSAP(x[1]-x[2]));
    out[2] = pSys[3]*SIGNAP(x[1]-x[2])*SQRT(f3*pSys[6]*ABSAP(x[1]-x[2]))+pSys[1]*x[4]-pSys[4]*(ctypeRNum)(4.8)*SQRT(f3*pSys[6]*x[2]);
    out[3] = u[0];
    out[4] = u[1];
    out[5] = u[2];
}
/** Jacobian df/dx multiplied by vector vec, i.e. (df/dx)^T*vec or vec^T*(df/dx) **/
void dfdx_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *vec, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
    ctypeRNum *pSys = (typeRNum*)userparam;
    ctypeRNum f1 = 1e4;
    ctypeRNum f2 = 1e-10;
    ctypeRNum f3 = 2;
    out[0] = vec[0]*((SQRT(f3))*f1*pSys[2]*x[5]*((POW2(TANH(f1*(x[0] - x[1])))) - 1)*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[0] - x[1]))))))))) - ((SQRT(f3))*pSys[2]*pSys[6]*x[5]*TANH(f1*(x[0] - x[1]))*(f3*x[0] - f3*x[1]))/ALMOSTZERO(4*(SQRT((f2 + (POW2((x[0] - x[1]))))))*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[0] - x[1]))))))))))) - vec[1]*((SQRT(f3))*f1*pSys[2]*x[5]*((POW2(TANH(f1*(x[0] - x[1])))) - 1)*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[0] - x[1]))))))))) - ((SQRT(f3))*pSys[2]*pSys[6]*x[5]*TANH(f1*(x[0] - x[1]))*(f3*x[0] - f3*x[1]))/ALMOSTZERO(4*(SQRT((f2 + (POW2((x[0] - x[1]))))))*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[0] - x[1])))))))))));

    out[1] = vec[1]*((SQRT(f3))*f1*pSys[3]*((POW2(TANH(f1*(x[1] - x[2])))) - 1)*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[1] - x[2]))))))))) + (SQRT(f3))*f1*pSys[2]*x[5]*((POW2(TANH(f1*(x[0] - x[1])))) - 1)*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[0] - x[1]))))))))) - ((SQRT(f3))*pSys[3]*pSys[6]*TANH(f1*(x[1] - x[2]))*(f3*x[1] - f3*x[2]))/ALMOSTZERO(4*(SQRT((f2 + (POW2((x[1] - x[2]))))))*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[1] - x[2])))))))))) - ((SQRT(f3))*pSys[2]*pSys[6]*x[5]*TANH(f1*(x[0] - x[1]))*(f3*x[0] - f3*x[1]))/ALMOSTZERO(4*(SQRT((f2 + (POW2((x[0] - x[1]))))))*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[0] - x[1]))))))))))) - vec[0]*((SQRT(f3))*f1*pSys[2]*x[5]*((POW2(TANH(f1*(x[0] - x[1])))) - 1)*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[0] - x[1]))))))))) - ((SQRT(f3))*pSys[2]*pSys[6]*x[5]*TANH(f1*(x[0] - x[1]))*(f3*x[0] - f3*x[1]))/ALMOSTZERO(4*(SQRT((f2 + (POW2((x[0] - x[1]))))))*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[0] - x[1]))))))))))) - vec[2]*((SQRT(f3))*f1*pSys[3]*((POW2(TANH(f1*(x[1] - x[2])))) - 1)*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[1] - x[2]))))))))) - ((SQRT(f3))*pSys[3]*pSys[6]*TANH(f1*(x[1] - x[2]))*(f3*x[1] - f3*x[2]))/ALMOSTZERO(4*(SQRT((f2 + (POW2((x[1] - x[2]))))))*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[1] - x[2])))))))))));

    out[2] = - vec[1]*((SQRT(f3))*f1*pSys[3]*((POW2(TANH(f1*(x[1] - x[2])))) - 1)*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[1] - x[2]))))))))) - ((SQRT(f3))*pSys[3]*pSys[6]*TANH(f1*(x[1] - x[2]))*(f3*x[1] - f3*x[2]))/ALMOSTZERO(4*(SQRT((f2 + (POW2((x[1] - x[2]))))))*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[1] - x[2]))))))))))) - vec[2]*((12*(SQRT(f3))*pSys[4]*pSys[6])/ALMOSTZERO(5*(SQRT((x[2]*pSys[6])))) - (SQRT(f3))*f1*pSys[3]*((POW2(TANH(f1*(x[1] - x[2])))) - 1)*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[1] - x[2]))))))))) + ((SQRT(f3))*pSys[3]*pSys[6]*TANH(f1*(x[1] - x[2]))*(f3*x[1] - f3*x[2]))/ALMOSTZERO(4*(SQRT((f2 + (POW2((x[1] - x[2]))))))*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[1] - x[2])))))))))));

    out[3] = pSys[0]*vec[0];

    out[4] = pSys[1]*vec[2];

    out[5] = (SQRT(f3))*pSys[2]*vec[1]*TANH(f1*(x[0] - x[1]))*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[0] - x[1]))))))))) - (SQRT(f3))*pSys[2]*vec[0]*TANH(f1*(x[0] - x[1]))*(SQRT((pSys[6]*(SQRT((f2 + (POW2((x[0] - x[1])))))))));
}
/** Jacobian df/du multiplied by vector vec, i.e. (df/du)^T*vec or vec^T*(df/du) **/
void dfdu_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *vec, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
    ctypeRNum *pSys = (ctypeRNum*)userparam;

    out[0] = vec[3];
    out[1] = vec[4];
    out[2] = vec[5];
}
/** Jacobian df/dp multiplied by vector vec, i.e. (df/dp)^T*vec or vec^T*(df/dp) **/
void dfdp_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *vec, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
}

/** Integral cost l(t,x(t),u(t),p,xdes,udes,userparam)
    -------------------------------------------------- **/
void lfct(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *xdes, ctypeRNum *udes, typeUSERPARAM *userparam)
{
    ctypeRNum* pSys = (ctypeRNum*)userparam;
    ctypeRNum* pCost = pSys + NPSYS; // Verschieben der Adresse in Zeiger p = userparam um NPSYS (=size(p)) nach vorn, um pCost zu erreichen wie in initData.m definiert

    out[0] = pCost[0] * POW2(x[0] - xdes[0])/(typeRNum)84100
        + pCost[1] * POW2(x[1] - xdes[1])/(typeRNum)84100
        + pCost[2] * POW2(x[2] - xdes[2])/(typeRNum)84100
        + pCost[6] * POW2(u[0] - udes[0])/(typeRNum)98
        + pCost[7] * POW2(u[1] - udes[1])/(typeRNum)98
        + pCost[8] * POW2(u[2] - udes[2])/(typeRNum)98;
}
/** Gradient dl/dx **/
void dldx(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *xdes, ctypeRNum *udes, typeUSERPARAM *userparam)
{
    ctypeRNum* pSys = (ctypeRNum*)userparam;
    ctypeRNum* pCost = pSys + NPSYS;

    out[0] = 2 * pCost[0] * (x[0] - xdes[0])/(typeRNum)84100;
    out[1] = 2 * pCost[1] * (x[1] - xdes[1])/(typeRNum)84100;
    out[2] = 2 * pCost[2] * (x[2] - xdes[2])/(typeRNum)84100;
    out[3] = 0;
    out[4] = 0;
    out[5] = 0;

}
/** Gradient dl/du **/
void dldu(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *xdes, ctypeRNum *udes, typeUSERPARAM *userparam)
{
    ctypeRNum* pSys = (ctypeRNum*)userparam;
    ctypeRNum* pCost = pSys + NPSYS;

    out[0] = 2 * pCost[6] * (u[0] - udes[0])/(typeRNum)98;
    out[1] = 2 * pCost[7] * (u[1] - udes[1])/(typeRNum)98;
    out[2] = 2 * pCost[8] * (u[2] - udes[2])/(typeRNum)98;
}
/** Gradient dl/dp **/
void dldp(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *xdes, ctypeRNum *udes, typeUSERPARAM *userparam)
{
}

/** Terminal cost V(T,x(T),p,xdes,userparam)
    ---------------------------------------- **/
void Vfct(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *xdes, typeUSERPARAM *userparam)
{
    ctypeRNum* pSys = (ctypeRNum*)userparam;
    ctypeRNum* pCost = pSys + NPSYS;

    out[0] = pCost[3] * POW2(x[0] - xdes[0])/(typeRNum)84100
        + pCost[4] * POW2(x[1] - xdes[1])/(typeRNum)84100
        + pCost[5] * POW2(x[2] - xdes[2])/(typeRNum)84100;

}
/** Gradient dV/dx **/
void dVdx(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *xdes, typeUSERPARAM *userparam)
{
    ctypeRNum* pSys = (ctypeRNum*)userparam;
    ctypeRNum* pCost = pSys + NPSYS;

    out[0] = 2 * pCost[3] * (x[0] - xdes[0])/(typeRNum)84100;
    out[1] = 2 * pCost[4] * (x[1] - xdes[1])/(typeRNum)84100;
    out[2] = 2 * pCost[5] * (x[2] - xdes[2])/(typeRNum)84100;
    out[3] = 0;
    out[4] = 0;
    out[5] = 0;
}
/** Gradient dV/dp **/
void dVdp(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *xdes, typeUSERPARAM *userparam)
{
}
/** Gradient dV/dT **/
void dVdT(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *xdes, typeUSERPARAM *userparam)
{
}

/** Equality constraints g(t,x(t),u(t),p,uperparam) = 0
    --------------------------------------------------- **/
void gfct(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Jacobian dg/dx multiplied by vector vec, i.e. (dg/dx)^T*vec or vec^T*(dg/dx) **/
void dgdx_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
/** Jacobian dg/du multiplied by vector vec, i.e. (dg/du)^T*vec or vec^T*(dg/du) **/
void dgdu_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
/** Jacobian dg/dp multiplied by vector vec, i.e. (dg/dp)^T*vec or vec^T*(dg/dp) **/
void dgdp_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}

/** Inequality constraints h(t,x(t),u(t),p,uperparam) <= 0
    ------------------------------------------------------ **/
void hfct(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
    ctypeRNum *pSys = (ctypeRNum*)userparam;
    ctypeRNum *pConstr = pSys + NPSYS + NPCOST;
    ctypeRNum *pOffset = pSys + NPSYS + NPCOST + NPCONSTR;
    ctypeRNum *pScale = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET;
    out[0] = ((x[0] - pOffset[0]) - pConstr[0])/pScale[0];
    out[1] = ((x[1] - pOffset[1]) - pConstr[1])/pScale[1];
    out[2] = ((x[2] - pOffset[2]) - pConstr[2])/pScale[2];
    out[3] = ((x[3] - pOffset[3]) - pConstr[3])/pScale[3];
    out[4] = ((x[4] - pOffset[4]) - pConstr[4])/pScale[4];
    out[5] = ((x[5] - pOffset[5]) - pConstr[5])/pScale[5];
    out[6] = (-(x[0] - pOffset[6]) + pConstr[6])/pScale[6];
    out[7] = (-(x[1] - pOffset[7]) + pConstr[7])/pScale[7];
    out[8] = (-(x[2] - pOffset[8]) + pConstr[8])/pScale[8];
    out[9] = (-(x[3] - pOffset[9]) + pConstr[9])/pScale[9];
    out[10] = (-(x[4] - pOffset[10]) + pConstr[10])/pScale[10];
    out[11] = (-(x[5] - pOffset[11]) + pConstr[11])/pScale[11];

    // Inequality for delta
    out[12] = ((u[0])-pConstr[12])/pScale[12];
    out[13] = ((u[1])-pConstr[13])/pScale[13];
    out[14] = ((u[2])-pConstr[14])/pScale[14];

    out[15] = (-(u[0])+pConstr[15])/pScale[15];
    out[16] = (-(u[1])+pConstr[16])/pScale[16];
    out[17] = (-(u[2])+pConstr[17])/pScale[17];

}
/** Jacobian dh/dx multiplied by vector vec, i.e. (dh/dx)^T*vec or vec^T*(dg/dx) **/
void dhdx_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
    ctypeRNum *pSys = (ctypeRNum*)userparam;
    ctypeRNum *pScale = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET;
    //[dhdx] = [nh x nx] = [12 x 6], tag_onenote:jacobi_dhdx
    out[0] = vec[0]/pScale[0]-vec[6]/pScale[6];
    out[1] = vec[1]/pScale[1]-vec[7]/pScale[7];
    out[2] = vec[2]/pScale[2]-vec[8]/pScale[8];
    out[3] = vec[3]/pScale[3]-vec[9]/pScale[9];
    out[4] = vec[4]/pScale[4]-vec[10]/pScale[10];
    out[5] = vec[5]/pScale[5]-vec[11]/pScale[11];

}
/** Jacobian dh/du multiplied by vector vec, i.e. (dh/du)^T*vec or vec^T*(dg/du) **/
void dhdu_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
    ctypeRNum *pSys = (ctypeRNum*)userparam;
    ctypeRNum *pScale = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET;
    out[0] = vec[12]/pScale[12]-vec[15]/pScale[15];
    out[1] = vec[13]/pScale[13]-vec[16]/pScale[16];
    out[2] = vec[14]/pScale[14]-vec[17]/pScale[17];

}
/** Jacobian dh/dp multiplied by vector vec, i.e. (dh/dp)^T*vec or vec^T*(dg/dp) **/
void dhdp_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}

/** Terminal equality constraints gT(T,x(T),p,uperparam) = 0
    -------------------------------------------------------- **/
void gTfct(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Jacobian dgT/dx multiplied by vector vec, i.e. (dgT/dx)^T*vec or vec^T*(dgT/dx) **/
void dgTdx_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
void tdgdu_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
/** Jacobian dgT/dp multiplied by vector vec, i.e. (dgT/dp)^T*vec or vec^T*(dgT/dp) **/
void dgTdp_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
/** Jacobian dgT/dT multiplied by vector vec, i.e. (dgT/dT)^T*vec or vec^T*(dgT/dT) **/
void dgTdT_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}

/** Terminal inequality constraints hT(T,x(T),p,uperparam) <= 0
    ----------------------------------------------------------- **/
void hTfct(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Jacobian dhT/dx multiplied by vector vec, i.e. (dhT/dx)^T*vec or vec^T*(dhT/dx) **/
void dhTdx_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
/** Jacobian dhT/dp multiplied by vector vec, i.e. (dhT/dp)^T*vec or vec^T*(dhT/dp) **/
void dhTdp_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
/** Jacobian dhT/dT multiplied by vector vec, i.e. (dhT/dT)^T*vec or vec^T*(dhT/dT) **/
void dhTdT_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}

/** Additional functions required for semi-implicit systems
    M*dx/dt(t) = f(t0+t,x(t),u(t),p) using the solver RODAS
    ------------------------------------------------------- **/
/** Jacobian df/dx in vector form (column-wise) **/
void dfdx(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Jacobian df/dx in vector form (column-wise) **/
void dfdxtrans(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Jacobian df/dt **/
void dfdt(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Jacobian d(dH/dx)/dt  **/
void dHdxdt(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *vec, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Mass matrix in vector form (column-wise, either banded or full matrix) **/
void Mfct(typeRNum *out, typeUSERPARAM *userparam)
{
}
/** Transposed mass matrix in vector form (column-wise, either banded or full matrix) **/
void Mtrans(typeRNum *out, typeUSERPARAM *userparam)
{
}

[Felix Mesmer]

Hi,

I downloaded your code and ran your example, but was not able to reproduce your results. For me the results looked good regardless of the prediction horizon (60, 90, 100). Are you using single or double precision floating points? Furthermore a lot less iterations were sufficient for me, i.e. 5/5 did perfectly fine.

Regarding your NaNs, in your defines you compare a floating point variable with zero. Here it would be better to make a comparison to a small epsilon. In addition, the macro SQRT() is already defined in "grampc_init.h", which could lead to problems.

I believe the error with the estimation of the penalty parameter is related to the NaNs you are seeing, since then the estimation result will also be NaN, which is not in the valid parameter range (0, inf).

Regards,
Felix

[hcl734]

Hi Felix,

thanks for your answer. I expressed myself unclearly, the simulation also runs fine for me with 60, 90 and 100. But for a horizon bigger than 112, e.g. 120 (N=240 in initData.m) the simulation terminates, see the attached zip.

I am using double precision with
grampc_init.h

#define USE_typeRNum   USE_DOUBLE /* USE_FLOAT */

I am using those defines as an approximation for discontinous functions like SIGN and ABS that are part of my problems ODEs. ALMOSTZERO I added to be sure not to do division by zero and the altered SQRT to avoid taking the root of a negative number.

#define POW2(a) ((a)*(a))
#define SIGNAP(x)                   tanh(10000*x)
#define ABSAP(x)                sqrt(pow(x,2)+1e-10)
#define TANH(x)                     tanh(x)
#define SQRT(x)           ((x) > 1e-15 ? sqrt(x): 1e-15)
#define ALMOSTZERO(x)           ((x) != 0 ? (x): 1e-16)

Is there a better way to approximate those in GRAMPC?

Best regards

[Felix Mesmer]

Hi,

I can reproduce your results/problem and I believe they are caused by numerical problems because of your system dynamic. To find out where exactly things go wrong you can attach Visual Studio to Matlab and debug your code.

One possible remedy could be to approximate the right hand side of your system around the critical points (e.g. x[0]-x[1] close to zero) with a polynomial, instead of approximating the sigmoid-function. By choosing the appropriate order you can ensure that your gradient is continuous.

Regards,
Felix

[hcl734]

Hi Felix,

thanks again for your advice.
I tried to debug with Visual Studio, building with startMPC([], 2, debug), attaching VS debugger to the process and then setting a break point in probfct_tank.c.
Unfortunately if I execute now with startMPC(1, 0) the breakpoint is not hit. Is probfct_tank.c not evaluated during simulation?
I validated that the debugging works for other (non GRAMPC) functions following these instructions:
https://de.mathworks.com/help/matlab/matlab_external/debugging-on-microsoft-windows-platforms.html
Is my workflow correct?

Appoximating the RHS around the critical points by polynomials you mean for example a taylor expansion around h1-h2 = 0 to use just when h1-h2 is close to zero and to use the "normal" dynamics when h1-h2 <<0 (h1-h2>>0) with the correct signs.
How would I implement something like this in the probfct.c?
Is a normal if-statement sufficient?

Best regrads

[Felix Mesmer]

Hi,

I encountered your problem with debugging too. You have to attach the VS debugger before the compiled mex-File is executed. The easiest way to achieve this is to start with "startMPC([],2,'debug')" and have a break point in Matlab in line 47 after the compilation.

A normal if-statement is sufficient, if you use it in the function for the RHS (i.e. out[0] = f_1(h1-h2)) and the corresponding gradients(i.e. the functions ffct, dfdx_vec, dfdu_vec, dfdp_vec). You can define a small epsilon, so that if |h1-h2| < epsilon the polynomial approximation f_p(h1-h2) is used. For the left and right border you have several equations to ensure that the transition from your function to the approximation is continuous (as well as for the corresponding partial derivatives), i.e. f_1(-epsilon) = f_p(-epsilon), f_1(+epsilon) = f_p(+epsilon) and the corresponding partial derivatives. With those equations you can now determine the coefficients of the polynomial.

Regards,
Felix

[hcl734]

Hi Felix,

great hint now I am able to debug with VS. I noticed that the first occurence of NaN-values is within the vec[i] variables when I change the xdes values for the setpoint change I am simulating. They are growing very fast and then become infinity or nan.
In which function are those values determined?

I took your advice and approximated the function around the critical points with a polynomial. However the problem with the NaN-values remains.
But I noticed that over the predicition the state values become negative (in violation with the constraints). Obviously my equations are not fit for negative values of the states.
So I introduced a direct constraint for the state values in probfct.c's ffct:

    if (x[0] < 1e-16) // Höhe im ersten Tank Null, damit kann Höhe im ersten Tank nur noch zunehmen nicht abnehmen
    {
        x[0] = 5e-17;
        if (out[0] < 1e-16)
        {
            out[0] = 1e-16;
        }
    }

so if x[0] is almost zero or smaller than zero it is set to a constant almost zero as well as the change of x[0] (out[0]) which can only be positive then. This represents the physical reality of my system where the state x[0] can't be smaller than 0.

However I am not sure how to set the dfdx_vec function accordingly.
If x goes to zero wouldn't dfdx go to infinity?

[Felix Mesmer]

Hi,
do you mean the vec values in the function “dfdx_vec”? In this case, they are the adjoint states and calculated in gampc_run.c in the function “evaluate_adjsys” (line 637).

Like Andreas wrote in the discussion “Constraints handling in Grampc”, the states in intermediate iterations can violate the constraints. Based on a gut feeling (based on you saying the predicted states become negative) I set user.opt.LineSearchMax = 1e-2, and then your simulation runs even with a prediction horizon of 150 s (Nhor = 300). If you set PLOT_STEPS = 1, PLOT_PAUSE = 1 in your startMPC and user.opt.MaxGradIter = 1 and user.opt.MaxMultIter = 1, you can see what happens without the change to the step size. The first step is too large and your tank heights become negative, which then leads to the NaNs in successive iterations.

Regards,
Felix

[hcl734]

Hi Felix,
setting LineSearchMax lower works great.
However it should be possible to account for the fact that negative values can occur in the problem formulation directly, no?
If I limit the states to positive values in ffct and dfdx before the actual equations with:

if (x[0] < 1e-16) 
    {
        x[0] = 5e-17;
    }

The states don't become negative anymore in the prediction plot however I see that they are still negative in grampc.rws.x after running the problem and the nan issue remains.
I attached a current version of my problem. Could you have a look at it?
I am not aware were the negative values in grampc.rws.x come from.

Best regards

[Felix Mesmer]

Hi,

it should definitely be possible, but you are not supposed to change the value of the state x from within the probfct.c. This is why all arguments besides out and userparam are constant...
What you have to do, is to change the values in out[0], if x[0] is close or equal to zero. I believe your state x[3] is some kind of in/outflow for you first state, i.e. if x[0] == 0, then x[3] < 0 is not possible and the respective part in out[0] = ... should be removed?

if ( (x[0] < 1e-16)  && (x[3] < 0.0) )
{
    out[0] = x[5]*pSys[2]*SIGNAP(x[0]-x[1])*SQRT2(f3*pSys[6]*ABSAP(x[0]-x[1]));
}
else
{
    out[0] = pSys[0]*x[3]-x[5]*pSys[2]*SIGNAP(x[0]-x[1])*SQRT2(f3*pSys[6]*ABSAP(x[0]-x[1]));
}

(The same is true for x[4] and out[2])

Regards,
Felix

[hcl734]

Hi Felix,

I see that it is problematic altering x[i] like this.
x[3] is indeed an inflow through a pump therefore it can only be positive so I would have to alter the condition to

if (  (x[3] < 0.0) )
{
    out[0] = x[5]*pSys[2]*SIGNAP(x[0]-x[1])*SQRT2(f3*pSys[6]*ABSAP(x[0]-x[1]));
}

and similarly I have to alter then to not allow negative heights in the tanks (x[0], x[1], x[2])

if (  (x[0] < 0.0) )
{
    out[0] = pSys[0]*x[3]-x[5]*pSys[2]*SIGNAP(x[1])*SQRT2(f3*pSys[6]*ABSAP(-x[1]));
}

and combine those conditions in their variations (e.g. x[3] < 0.0 && x[0] < 0.0).
I summarized this into defining some floats for the states and setting them zero for negative values.
And added the condition that x[4], x[5], x[6] can never be negative.

    double X0, X1, X2, X3, X4, X5, U0, U1, U2;

    // C arrays werden ab 0 indiziert --> p(1) = pSys[0]...; p(7) = pSys[6]; p(6) = pSys[5]
    X0 = x[0];
    X1 = x[1];
    X2 = x[2];
    X3 = x[3];
    X4 = x[4];
    X5 = x[5];
    U0 = u[0];
    U1 = u[1];
    U2 = u[2];
    if (( x[0] < 1 ))
    {
       X0 = 0.0; 
    }
    if (( x[1] < 1 ))
    {
       X1 = 0.0; 
    }
    if (( x[2] < 1 ))
    {
       X2 = 0.0; 
    }
    if (( x[3] < 0.1 ))
    {
       X3 = 0.0;
       U0 = 0.0;
    }
    if (( x[4] < 0.1 ))
    {
       X4 = 0.0; 
       U1= 0.0;

    }
    if (( x[5] < 0.1 ))
    {
       X5 = 0.0; 
       U2= 0.0;
    }

    out[0] = pSys[0]*X3-X5*pSys[2]*SIGNAP(X0-X1)*SQRT(f3*pSys[6]*ABSAP(X0-X1));
    out[1] = X5*pSys[2]*SIGNAP(X0-X1)*SQRT(f3*pSys[6]*ABSAP(X0-X1))-pSys[3]*SIGNAP(X1-X2)*SQRT(f3*pSys[6]*ABSAP(X1-X2));
    out[2] = pSys[3]*SIGNAP(X1-X2)*SQRT(f3*pSys[6]*ABSAP(X1-X2))+pSys[1]*X4-pSys[4]*(ctypeRNum)(4.8)*SQRT(f3*pSys[6]*X2);

    out[3] = U0;
    out[4] = U1;
    out[5] = U2;

And indeed now the solution in grampc.rws.x is always positive.
But I still get the NaN error.
What am I doing wrong here?

[Felix Mesmer]

If you do it like this, your system dynamic will not be continuous. This might lead to numerical problems in the integration of your adjoint states. You didn’t write so explicitly, but I assume you also added this change to every function in your probfct-file?
In addition, your predicted states 4-6 can very well be negative, since

    out[3] = U0;
    out[4] = U1;
    out[5] = U2;

Like I wrote earlier, after changing the option LineSearchMax, I didn’t have any problems running your example (no NaNs appearing), did this not work well enough for you?

[hcl734]

Just to close this thread. There was a mistake in my macro definition:

#define SIGNAP(x)                   tanh(10000*x)
 // --> SIGNAP(x+y) = tanh(10000*x+ y

which caused the numerical issues.

I experimented with limiting the states with sigmoid functions to avoid the non-continuity that you mentioned which of course was computationally expensive. In the end I limited the functions with the if clauses as described before and of course also adding a restriction on the change of U0 as shown below.
Now I do not have problems with negative values anymore and the results for the experiment are good.
Thanks for your help!

problfct.c

// Formulierung mit der Approximation durch tanh und harte Begrenzungen

#include "probfct.h"

#define NPSYS 7 // size(p), s.u.
#define NPCOST 9 // size(p), s.u.
#define NPCONSTR 18
#define NPOFFSET 18
#define NPSCALE 18
#define NPTS 1

#define cAP 1e-5
 /* square macro */
#define POW2(a) ((a)*(a))
#define SQRT(x)           ((x) > 1e-15 ? sqrt(x): 1e-16)

#define POW2(a) ((a)*(a))
#define POW3(a) ((a)*(a)*(a))
#define POW4(a) ((a)*(a)*(a)*(a))
#define POW5(a) ((a)*(a)*(a)*(a)*(a))
#define POW7(a) ((a)*(a)*(a)*(a)*(a)*(a)*(a))
#define ex32(x)                 POW3(SQRT(x))
#define ex12(x)     ((x) > 1e-15 ? sqrt(x): 1e-15)
#define ex42(x)     POW4(SQRT(x))
#define ex52(x)     POW5(SQRT(x))
#define ex72(x)     POW7(SQRT(x))

/** OCP dimensions: states (Nx), controls (Nu), parameters (Np), equalities (Ng),
    inequalities (Nh), terminal equalities (NgT), terminal inequalities (NhT) **/
void ocp_dim(typeInt *Nx, typeInt *Nu, typeInt *Np, typeInt *Ng, typeInt *Nh, typeInt *NgT, typeInt *NhT, typeUSERPARAM *userparam)
{
    *Nx = 6; // h1=X0 h2=X1 h3=X2 p1=X3 p2=X4 v12=X5
    *Nu = 3; // delta_p1=U0 delta_p2=U1 delta_v12=U2
    *Np = 0; //
    *Nh = 18; //
    *Ng = 0;
    *NgT = 0;
    *NhT = 0;
}

/** System function f(t,x,u,p,userparam)
    ------------------------------------ **/
void ffct(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
    // typeRNum ist double
    // ctypeRNum ist constant double
    ctypeRNum *pSys = (typeRNum*)userparam;
    ctypeRNum f3 = 2;
    ctypeRNum *Ts = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET + NPSCALE;
    ctypeRNum *dev = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET + NPSCALE + NPTS;

    double X0, X1, X2, X3, X4, X5, U0, U1, U2, U2cal, U2new, delta_h1, delta_h2, delta_h3, delta_h1_func, delta_h2_func, delta_h3_func, test1, test2, test3, test4;
    X0 = x[0];
    X1 = x[1];
    X2 = x[2];
    X3 = x[3];
    X4 = x[4];
    X5 = x[5];
    U0 = u[0];
    U1 = u[1];
    U2 = u[2];

    if (x[0] < 1e-16) 
    {
        X0 = 5e-17;
    }
    else if(x[0] > 300) 
    {
        X0 = 300;
    }

    if (x[1] < 1e-16) 
    {
        X1 = 5e-17;
    }
    else if (x[1] > 300) 
    {
        X1 = 300;
    }

    if (x[2] < 1e-16) 
    {
        X2 = 5e-17;
    }
    else if (x[2] > 300) 
    {
        X2 = 300;
    }

    if (x[3] < 1e-16) 
    {
        X3 = 5e-17;
    }
    else if (x[3] > 9.9) 
    {
        X3 = 9.9;
    }

    if (x[4] < 1e-16) 
    {
        X4 = 5e-17;
    }
    else if (x[4] > 9.9) 
    {
        X4 = 9.9;
    }

    if (x[5] < 1e-16) 
    {
        X5 = 5e-17;
    }
    else if (x[5] > 9.9) 
    {
        X5 = 9.9;
    }    

    if (U0*Ts[0]+x[3] > 9.9) 
    {
        U0 = -5e-17;
    }
    else if (U0*Ts[0]+x[3] < 0) 
    {
        U0 = 5e-17;
    }

    if (U1*Ts[0]+x[4] > 9.9) 
    {
        U1 = -5e-17;
    }
    else if (U1*Ts[0]+x[4] < 0) 
    {
        U1 = 5e-17;
    }

    if (U2*Ts[0]+x[5] > 9.9) 
    {
        U2 = -5e-17;
    }
    else if (U2*Ts[0]+x[5] < 0) 
    {
        U2 = 5e-17;
    }

    delta_h1 = X3*pSys[0] - ((ex12(2))*pSys[2]*X5*(X0 - X1)*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1))))))))))/(ex12((cAP + (POW2((X0 - X1))))))+ dev[0];

    delta_h2 = ((ex12(2))*pSys[2]*X5*(X0 - X1)*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1))))))))))/(ex12((cAP + (POW2((X0 - X1)))))) - ((ex12(2))*pSys[3]*(X1 - X2)*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2))))))))))/(ex12((cAP + (POW2((X1 - X2))))))+ dev[1];

    delta_h3 = X4*pSys[1] - (24*(ex12(2))*pSys[4]*(ex12((X2*pSys[6]))))/5 + ((ex12(2))*pSys[3]*(X1 - X2)*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2))))))))))/(ex12((cAP + (POW2((X1 - X2))))))+ dev[2];

    if (delta_h1*Ts[0]+x[0] > 300) 
    {
        delta_h1= -5e-17;
    }
    else if (delta_h1*Ts[0]+x[0] < 0) 
    {
        delta_h1 = 5e-17;
    }

    if (delta_h2*Ts[0]+x[1] > 300) 
    {
        delta_h2 = -5e-17;
    }
    else if (delta_h2*Ts[0]+x[1] < 0) 
    {
        delta_h2 = 5e-17;
    }

    if (delta_h3*Ts[0]+x[2] > 300) 
    {
        delta_h3 = -5e-17;
    }
    else if (delta_h3*Ts[0]+x[2] < 0) 
    {
        delta_h3 = 5e-17;
    }

    out[0] = delta_h1;
    out[1] = delta_h2;
    out[2] = delta_h3;
    out[3] = U0;
    out[4] = U1;
    out[5] = U2;
}
/** Jacobian df/dx multiplied by vector vec, i.e. (df/dx)^T*vec or vec^T*(df/dx) **/
void dfdx_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *vec, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
    ctypeRNum *pSys = (typeRNum*)userparam;
    ctypeRNum f1 = 1e4;
    ctypeRNum f2 = 1e-10;
    ctypeRNum f3 = 2;
    ctypeRNum *Ts = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET + NPSCALE;

    double X0, X1, X2, X3, X4, X5, U0, U1, U2, U2cal, U2new;
    X0 = x[0];
    X1 = x[1];
    X2 = x[2];
    X3 = x[3];
    X4 = x[4];
    X5 = x[5];
    U0 = u[0];
    U1 = u[1];
    U2 = u[2];

    if (x[0] < 1e-16) 
    {
        X0 = 5e-17;
    }
    else if(x[0] > 300) 
    {
        X0 = 300;
    }

    if (x[1] < 1e-16) 
    {
        X1 = 5e-17;
    }
    else if (x[1] > 300) 
    {
        X1 = 300;
    }

    if (x[2] < 1e-16) 
    {
        X2 = 5e-17;
    }
    else if (x[2] > 300) 
    {
        X2 = 300;
    }

    if (x[3] < 1e-16) 
    {
        X3 = 5e-17;
    }
    else if (x[3] > 9.9) 
    {
        X3 = 9.9;
    }

    if (x[4] < 1e-16) 
    {
        X4 = 5e-17;
    }
    else if (x[4] > 9.9) 
    {
        X4 = 9.9;
    }

    if (x[5] < 1e-16) 
    {
        X5 = 5e-17;
    }
    else if (x[5] > 9.9) 
    {
        X5 = 9.9;
    }    

    if (U0*Ts[0]+x[3] > 9.9) 
    {
        U0 = -5e-17;
    }
    else if (U0*Ts[0]+x[3] < 0) 
    {
        U0 = 5e-17;
    }

    if (U1*Ts[0]+x[4] > 9.9) 
    {
        U1 = -5e-17;
    }
    else if (U1*Ts[0]+x[4] < 0) 
    {
        U1 = 5e-17;
    }

    if (U2*Ts[0]+x[5] > 9.9) 
    {
        U2 = -5e-17;
    }
    else if (U2*Ts[0]+x[5] < 0) 
    {
        U2 = 5e-17;
    }
out[0] = vec[1]*(((ex12(2))*pSys[2]*X5*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1))))))))))/(ex12((cAP + (POW2((X0 - X1)))))) - ((ex12(2))*pSys[2]*X5*(X0 - X1)*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1)))))))))*(2*X0 - 2*X1))/(2*(ex32((cAP + (POW2((X0 - X1))))))) + ((ex12(2))*pSys[2]*pSys[6]*X5*(X0 - X1)*(2*X0 - 2*X1))/(4*(cAP + (POW2((X0 - X1))))*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1))))))))))) - vec[0]*(((ex12(2))*pSys[2]*X5*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1))))))))))/(ex12((cAP + (POW2((X0 - X1)))))) - ((ex12(2))*pSys[2]*X5*(X0 - X1)*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1)))))))))*(2*X0 - 2*X1))/(2*(ex32((cAP + (POW2((X0 - X1))))))) + ((ex12(2))*pSys[2]*pSys[6]*X5*(X0 - X1)*(2*X0 - 2*X1))/(4*(cAP + (POW2((X0 - X1))))*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1)))))))))));
out[1] = vec[2]*(((ex12(2))*pSys[3]*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2))))))))))/(ex12((cAP + (POW2((X1 - X2)))))) - ((ex12(2))*pSys[3]*(X1 - X2)*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2)))))))))*(2*X1 - 2*X2))/(2*(ex32((cAP + (POW2((X1 - X2))))))) + ((ex12(2))*pSys[3]*pSys[6]*(X1 - X2)*(2*X1 - 2*X2))/(4*(cAP + (POW2((X1 - X2))))*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2))))))))))) - vec[1]*(((ex12(2))*pSys[3]*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2))))))))))/(ex12((cAP + (POW2((X1 - X2)))))) + ((ex12(2))*pSys[2]*X5*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1))))))))))/(ex12((cAP + (POW2((X0 - X1)))))) - ((ex12(2))*pSys[3]*(X1 - X2)*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2)))))))))*(2*X1 - 2*X2))/(2*(ex32((cAP + (POW2((X1 - X2))))))) + ((ex12(2))*pSys[3]*pSys[6]*(X1 - X2)*(2*X1 - 2*X2))/(4*(cAP + (POW2((X1 - X2))))*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2)))))))))) - ((ex12(2))*pSys[2]*X5*(X0 - X1)*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1)))))))))*(2*X0 - 2*X1))/(2*(ex32((cAP + (POW2((X0 - X1))))))) + ((ex12(2))*pSys[2]*pSys[6]*X5*(X0 - X1)*(2*X0 - 2*X1))/(4*(cAP + (POW2((X0 - X1))))*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1))))))))))) + vec[0]*(((ex12(2))*pSys[2]*X5*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1))))))))))/(ex12((cAP + (POW2((X0 - X1)))))) - ((ex12(2))*pSys[2]*X5*(X0 - X1)*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1)))))))))*(2*X0 - 2*X1))/(2*(ex32((cAP + (POW2((X0 - X1))))))) + ((ex12(2))*pSys[2]*pSys[6]*X5*(X0 - X1)*(2*X0 - 2*X1))/(4*(cAP + (POW2((X0 - X1))))*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1)))))))))));
out[2] = vec[1]*(((ex12(2))*pSys[3]*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2))))))))))/(ex12((cAP + (POW2((X1 - X2)))))) - ((ex12(2))*pSys[3]*(X1 - X2)*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2)))))))))*(2*X1 - 2*X2))/(2*(ex32((cAP + (POW2((X1 - X2))))))) + ((ex12(2))*pSys[3]*pSys[6]*(X1 - X2)*(2*X1 - 2*X2))/(4*(cAP + (POW2((X1 - X2))))*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2))))))))))) - vec[2]*((12*(ex12(2))*pSys[4]*pSys[6])/(5*(ex12((X2*pSys[6])))) + ((ex12(2))*pSys[3]*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2))))))))))/(ex12((cAP + (POW2((X1 - X2)))))) - ((ex12(2))*pSys[3]*(X1 - X2)*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2)))))))))*(2*X1 - 2*X2))/(2*(ex32((cAP + (POW2((X1 - X2))))))) + ((ex12(2))*pSys[3]*pSys[6]*(X1 - X2)*(2*X1 - 2*X2))/(4*(cAP + (POW2((X1 - X2))))*(ex12((pSys[6]*(ex12((cAP + (POW2((X1 - X2)))))))))));
out[3] = pSys[0]*vec[0];
out[4] = pSys[1]*vec[2];
out[5] = ((ex12(2))*pSys[2]*vec[1]*(X0 - X1)*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1))))))))))/(ex12((cAP + (POW2((X0 - X1)))))) - ((ex12(2))*pSys[2]*vec[0]*(X0 - X1)*(ex12((pSys[6]*(ex12((cAP + (POW2((X0 - X1))))))))))/(ex12((cAP + (POW2((X0 - X1))))));
}
/** Jacobian df/du multiplied by vector vec, i.e. (df/du)^T*vec or vec^T*(df/du) **/
void dfdu_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *vec, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
    ctypeRNum *pSys = (ctypeRNum*)userparam;
    ctypeRNum *Ts = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET + NPSCALE;

    double X0, X1, X2, X3, X4, X5, U0, U1, U2, U2cal, U2new;
    X0 = x[0];
    X1 = x[1];
    X2 = x[2];
    X3 = x[3];
    X4 = x[4];
    X5 = x[5];
    U0 = u[0];
    U1 = u[1];
    U2 = u[2];

out[0] = vec[3];
out[1] = vec[4];
out[2] = vec[5];

}
/** Jacobian df/dp multiplied by vector vec, i.e. (df/dp)^T*vec or vec^T*(df/dp) **/
void dfdp_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *vec, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
}

/** Integral cost l(t,x(t),u(t),p,xdes,udes,userparam)
    -------------------------------------------------- **/
void lfct(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *xdes, ctypeRNum *udes, typeUSERPARAM *userparam)
{
    ctypeRNum* pSys = (ctypeRNum*)userparam;
    ctypeRNum* pCost = pSys + NPSYS; // Verschieben der Adresse in Zeiger p = userparam um NPSYS (=size(p)) nach vorn, um pCost zu erreichen wie in initData.m definiert
    ctypeRNum *Ts = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET + NPSCALE;

    double X0, X1, X2, X3, X4, X5, U0, U1, U2, U2cal;
    X0 = x[0];
    X1 = x[1];
    X2 = x[2];
    X3 = x[3];
    X4 = x[4];
    X5 = x[5];
    U0 = u[0];
    U1 = u[1];
    U2 = u[2];

    out[0] = pCost[0] * POW2(X0 - xdes[0])//(typeRNum)84100
        + pCost[1] * POW2(X1 - xdes[1])//(typeRNum)84100
        + pCost[2] * POW2(X2 - xdes[2])//(typeRNum)84100
        + pCost[6] * POW2(U0 - udes[0])//(typeRNum)98
        + pCost[7] * POW2(U1 - udes[1])//(typeRNum)98
        + pCost[8] * POW2(U2 - udes[2]);//(typeRNum)98;
}
/** Gradient dl/dx **/
void dldx(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *xdes, ctypeRNum *udes, typeUSERPARAM *userparam)
{
    ctypeRNum* pSys = (ctypeRNum*)userparam;
    ctypeRNum* pCost = pSys + NPSYS;
    ctypeRNum *Ts = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET + NPSCALE;

    double X0, X1, X2, X3, X4, X5, U0, U1, U2, U2cal;
    X0 = x[0];
    X1 = x[1];
    X2 = x[2];
    X3 = x[3];
    X4 = x[4];
    X5 = x[5];
    U0 = u[0];
    U1 = u[1];
    U2 = u[2];

    out[0] = 2 * pCost[0] * (X0 - xdes[0]);//(typeRNum)84100;
    out[1] = 2 * pCost[1] * (X1 - xdes[1]);//(typeRNum)84100;
    out[2] = 2 * pCost[2] * (X2 - xdes[2]);//(typeRNum)84100;
    out[3] = 0;
    out[4] = 0;
    out[5] = 0;

}
/** Gradient dl/du **/
void dldu(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *xdes, ctypeRNum *udes, typeUSERPARAM *userparam)
{
    ctypeRNum* pSys = (ctypeRNum*)userparam;
    ctypeRNum* pCost = pSys + NPSYS;
    ctypeRNum *Ts = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET + NPSCALE;

    double X0, X1, X2, X3, X4, X5, U0, U1, U2, U2cal;
    X0 = x[0];
    X1 = x[1];
    X2 = x[2];
    X3 = x[3];
    X4 = x[4];
    X5 = x[5];
    U0 = u[0];
    U1 = u[1];
    U2 = u[2];

    out[0] = 2 * pCost[6] * (U0 - udes[0]);//(typeRNum)98;
    out[1] = 2 * pCost[7] * (U1 - udes[1]);//(typeRNum)98;
    out[2] = 2 * pCost[8] * (U2 - udes[2]);//(typeRNum)98;
}
/** Gradient dl/dp **/
void dldp(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *xdes, ctypeRNum *udes, typeUSERPARAM *userparam)
{
}

/** Terminal cost V(T,x(T),p,xdes,userparam)
    ---------------------------------------- **/
void Vfct(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *xdes, typeUSERPARAM *userparam)
{
    ctypeRNum* pSys = (ctypeRNum*)userparam;
    ctypeRNum* pCost = pSys + NPSYS;

    double X0, X1, X2, X3, X4, X5;
    X0 = x[0];
    X1 = x[1];
    X2 = x[2];
    X3 = x[3];
    X4 = x[4];
    X5 = x[5];

    out[0] = pCost[3] * POW2(X0 - xdes[0])//(typeRNum)84100
        + pCost[4] * POW2(X1 - xdes[1])//(typeRNum)84100
        + pCost[5] * POW2(X2 - xdes[2]);//(typeRNum)84100;

}
/** Gradient dV/dx **/
void dVdx(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *xdes, typeUSERPARAM *userparam)
{
    ctypeRNum* pSys = (ctypeRNum*)userparam;
    ctypeRNum* pCost = pSys + NPSYS;

    double X0, X1, X2, X3, X4, X5;
    X0 = x[0];
    X1 = x[1];
    X2 = x[2];
    X3 = x[3];
    X4 = x[4];
    X5 = x[5];

    out[0] = 2 * pCost[3] * (X0 - xdes[0]);//(typeRNum)84100;
    out[1] = 2 * pCost[4] * (X1 - xdes[1]);//(typeRNum)84100;
    out[2] = 2 * pCost[5] * (X2 - xdes[2]);//(typeRNum)84100;
    out[3] = 0;
    out[4] = 0;
    out[5] = 0;
}
/** Gradient dV/dp **/
void dVdp(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *xdes, typeUSERPARAM *userparam)
{
}
/** Gradient dV/dT **/
void dVdT(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *xdes, typeUSERPARAM *userparam)
{
}

/** Equality constrainTs[0] g(t,x(t),u(t),p,uperparam) = 0
    --------------------------------------------------- **/
void gfct(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Jacobian dg/dx multiplied by vector vec, i.e. (dg/dx)^T*vec or vec^T*(dg/dx) **/
void dgdx_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
/** Jacobian dg/du multiplied by vector vec, i.e. (dg/du)^T*vec or vec^T*(dg/du) **/
void dgdu_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
/** Jacobian dg/dp multiplied by vector vec, i.e. (dg/dp)^T*vec or vec^T*(dg/dp) **/
void dgdp_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}

/** Inequality constraints h(t,x(t),u(t),p,uperparam) <= 0
    ------------------------------------------------------ **/
void hfct(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
    ctypeRNum *pSys = (ctypeRNum*)userparam;
    ctypeRNum *pConstr = pSys + NPSYS + NPCOST;
    ctypeRNum *pOffset = pSys + NPSYS + NPCOST + NPCONSTR;
    ctypeRNum *pScale = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET;
    ctypeRNum *Ts = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET + NPSCALE;

    double X0, X1, X2, X3, X4, X5, U0, U1, U2, U2cal;
    X0 = x[0];
    X1 = x[1];
    X2 = x[2];
    X3 = x[3];
    X4 = x[4];
    X5 = x[5];
    U0 = u[0];
    U1 = u[1];
    U2 = u[2];

    out[0] = ((X0 - pOffset[0]) - pConstr[0])/pScale[0];
    out[1] = ((X1 - pOffset[1]) - pConstr[1])/pScale[1];
    out[2] = ((X2 - pOffset[2]) - pConstr[2])/pScale[2];
    out[3] = ((X3 - pOffset[3]) - pConstr[3])/pScale[3];
    out[4] = ((X4 - pOffset[4]) - pConstr[4])/pScale[4];
    out[5] = ((X5 - pOffset[5]) - pConstr[5])/pScale[5];
    out[6] = (-(X0 - pOffset[6]) + pConstr[6])/pScale[6];
    out[7] = (-(X1 - pOffset[7]) + pConstr[7])/pScale[7];
    out[8] = (-(X2 - pOffset[8]) + pConstr[8])/pScale[8];
    out[9] = (-(X3 - pOffset[9]) + pConstr[9])/pScale[9];
    out[10] = (-(X4 - pOffset[10]) + pConstr[10])/pScale[10];
    out[11] = (-(X5 - pOffset[11]) + pConstr[11])/pScale[11];

    // Inequality for delta
    out[12] = ((U0)-pConstr[12])/pScale[12];
    out[13] = ((U1)-pConstr[13])/pScale[13];
    out[14] = ((U2)-pConstr[14])/pScale[14];

    out[15] = (-(U0)+pConstr[15])/pScale[15];
    out[16] = (-(U1)+pConstr[16])/pScale[16];
    out[17] = (-(U2)+pConstr[17])/pScale[17];

}
/** Jacobian dh/dx multiplied by vector vec, i.e. (dh/dx)^T*vec or vec^T*(dg/dx) **/
void dhdx_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
    ctypeRNum *pSys = (ctypeRNum*)userparam;
    ctypeRNum *pScale = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET;
    //[dhdx] = [nh x nx] = [12 x 6], tag_onenote:jacobi_dhdx
    ctypeRNum *Ts = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET + NPSCALE;

    out[0] = vec[0]/pScale[0]-vec[6]/pScale[6];
    out[1] = vec[1]/pScale[1]-vec[7]/pScale[7];
    out[2] = vec[2]/pScale[2]-vec[8]/pScale[8];
    out[3] = vec[3]/pScale[3]-vec[9]/pScale[9];
    out[4] = vec[4]/pScale[4]-vec[10]/pScale[10];
    out[5] = vec[5]/pScale[5]-vec[11]/pScale[11];

}
/** Jacobian dh/du multiplied by vector vec, i.e. (dh/du)^T*vec or vec^T*(dg/du) **/
void dhdu_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{

    ctypeRNum *pSys = (ctypeRNum*)userparam;
    ctypeRNum *pScale = pSys + NPSYS + NPCOST + NPCONSTR + NPOFFSET;
    out[0] = vec[12]/pScale[12]-vec[15]/pScale[15];
    out[1] = vec[13]/pScale[13]-vec[16]/pScale[16];
    out[2] = vec[14]/pScale[14]-vec[17]/pScale[17];

}
/** Jacobian dh/dp multiplied by vector vec, i.e. (dh/dp)^T*vec or vec^T*(dg/dp) **/
void dhdp_vec(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}

/** Terminal equality constraints gT(T,x(T),p,uperparam) = 0
    -------------------------------------------------------- **/
void gTfct(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Jacobian dgT/dx multiplied by vector vec, i.e. (dgT/dx)^T*vec or vec^T*(dgT/dx) **/
void dgTdx_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
void tdgdu_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
/** Jacobian dgT/dp multiplied by vector vec, i.e. (dgT/dp)^T*vec or vec^T*(dgT/dp) **/
void dgTdp_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
/** Jacobian dgT/dT multiplied by vector vec, i.e. (dgT/dT)^T*vec or vec^T*(dgT/dT) **/
void dgTdT_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}

/** Terminal inequality constraints hT(T,x(T),p,uperparam) <= 0
    ----------------------------------------------------------- **/
void hTfct(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Jacobian dhT/dx multiplied by vector vec, i.e. (dhT/dx)^T*vec or vec^T*(dhT/dx) **/
void dhTdx_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
/** Jacobian dhT/dp multiplied by vector vec, i.e. (dhT/dp)^T*vec or vec^T*(dhT/dp) **/
void dhTdp_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}
/** Jacobian dhT/dT multiplied by vector vec, i.e. (dhT/dT)^T*vec or vec^T*(dhT/dT) **/
void dhTdT_vec(typeRNum *out, ctypeRNum T, ctypeRNum *x, ctypeRNum *p, ctypeRNum *vec, typeUSERPARAM *userparam)
{
}

/** Additional functions required for semi-implicit systems
    M*dx/dt(t) = f(t0+t,x(t),u(t),p) using the solver RODAS
    ------------------------------------------------------- **/
/** Jacobian df/dx in vector form (column-wise) **/
void dfdx(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Jacobian df/dx in vector form (column-wise) **/
void dfdxtrans(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Jacobian df/dt **/
void dfdt(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Jacobian d(dH/dx)/dt  **/
void dHdxdt(typeRNum *out, ctypeRNum t, ctypeRNum *x, ctypeRNum *u, ctypeRNum *vec, ctypeRNum *p, typeUSERPARAM *userparam)
{
}
/** Mass matrix in vector form (column-wise, either banded or full matrix) **/
void Mfct(typeRNum *out, typeUSERPARAM *userparam)
{
}
/** Transposed mass matrix in vector form (column-wise, either banded or full matrix) **/
void Mtrans(typeRNum *out, typeUSERPARAM *userparam)
{
}