function [x,f,eflag,outpt,p0] = runobjconst(x0,lb,ub,opts,p1,p2,p3,fobj)
    % This function runs the optimization problem and sets it up for a more
    % efficient use of parallel computation. 
    % Inputs:
    %   x0: Initial point
    %   lb: Lower bound
    %   ub: Upper bound
    %   opts: Algorithm options
    %   Perturbations: 
    %       p1: Increment in the inlet temperature of BOG [K]
    %       p2: Modifier for the flow of BOG
    %       p3: Modifier for the incoming environment heat
    %   fobj: optimization mode:
    %       fobj = 1: normal optimzation
    %       fobj = 2: dummy objective function (simulation)
    %
    % Outputs:
    %   x: Optimal point
    %   f: Optimal objective function value
    %   eflag: Optimizers flag
    %   outpt: Optimizer output summary
    %   p0: point at iteration nr 10 for future initialization
    %% Initial parameter definition
    
    if nargin == 1
        opts = [];
        lb = [];
        ub = [];
    end
    xLast = [];
    myf = []; % Use for objective at xLast
    myc = []; % Use for nonlinear inequality constraint
    myceq = []; % Use for nonlinear equality constraint
    p0 = []; % Intermediate state for next iterations
    opts.OutputFcn=@outfun;
    
    
    fun = @objfun; % the objective function, nested below
    cfun = @constr; % the constraint function, nested below
    
    %% Optimization
    switch fobj
        case 1
            [x,f,eflag,outpt] = fmincon( @(w)fun(w,p1,p2,p3)  ,x0,[],[],[],[],lb,ub,@(w)cfun(w,p1,p2,p3),opts);
        case 2
            [x,f,eflag,outpt] = fmincon( '1',x0,[],[],[],[],lb,ub,@(w)cfun(w,p1,p2,p3),opts);
    end    
    
    %% Definition of the objective function based on the model
    function y = objfun(x,p1,p2,p3)
        if ~isequal(x,xLast) % Check if computation is necessary
            [myf,myc,myceq] = MainEO(x,p1,p2,p3);
            xLast = x;
        end
        y = myf;
    end
    %% Definition of the constraints based on the model
    function [c, ceq] = constr(x,p1,p2,p3)
        if ~isequal(x,xLast) % Check if computation is necessary
            [myf,myc,myceq] = MainEO(x,p1,p2,p3);
            xLast = x;
        end
        c = myc;
        ceq = myceq;
    end
    %% Function to extract an intermediate point data
    function stop = outfun(x,optimValues,state)
        stop=false;
        if optimValues.iteration == 15
            p0 = x;
        end
    end
end