new file mode 100644
--- /dev/null
+++ b/scripts/optimization/glpk.m
@@ -0,0 +1,309 @@
+function varargout=glpk(varargin)
+## GLPK - An Octave Interface for the GNU GLPK library
+##
+## This routine calls the glpk library to solve an LP/MIP problem. A typical
+## LP problem has following structure:
+##
+##           [min|max] C'x
+##            s.t.
+##                 Ax ["="|"<="|">="] b
+##                 {x <= UB}
+##                 {x >= LB}
+##
+## The calling syntax is:
+## [XOPT,FOPT,STATUS,EXTRA]=glpkmex(SENSE,C,...
+##                                      A,B,CTYPE,LB,UB,...
+##                                      VARTYPE,PARAM,LPSOLVER,SAVE)
+##
+## For a quick reference to the syntax just type glpk at command prompt.
+##
+## The minimum number of input arguments is 4 (SENSE,C,A,B). In this case we
+## assume all the constraints are '<=' and all the variables are continuous.
+##
+## --- INPUTS ---
+##
+## SENSE:     indicates whether the problem is a minimization
+##            or maximization problem.
+##            SENSE = 1 minimize
+##            SENSE = -1 maximize.
+##
+## C:         A column array containing the objective function
+##            coefficients.
+##
+## A:         A matrix containing the constraints coefficients.
+##
+## B:         A column array containing the right-hand side value for
+##            each constraint in the constraint matrix.
+##
+## CTYPE:     A column array containing the sense of each constraint
+##            in the constraint matrix.
+##            CTYPE(i) = 'F'  Free (unbounded) variable
+##            CTYPE(i) = 'U'  "<=" Variable with upper bound
+##            CTYPE(i) = 'S'  "="  Fixed Variable
+##            CTYPE(i) = 'L'  ">=" Variable with lower bound
+##            CTYPE(i) = 'D'  Double-bounded variable
+##            (This is case sensitive).
+##
+## LB:        An array of at least length numcols containing the lower
+##            bound on each of the variables.
+##
+## UB:        An array of at least length numcols containing the upper
+##            bound on each of the variables.
+##
+## VARTYPE:   A column array containing the types of the variables.
+##            VARTYPE(i) = 'C' continuous variable
+##            VARTYPE(i) = 'I' Integer variable
+##            (This is case sensitive).
+##
+## PARAM:     A structure containing some parameters used to define
+##            the behavior of solver. For more details type
+##            HELP GLPKPARAMS.
+##
+## LPSOLVER:  Selects which solver using to solve LP problems.
+##            LPSOLVER=1  Revised Simplex Method
+##            LPSOLVER=2  Interior Point Method
+##            If the problem is a MIP problem this flag will be ignored.
+##
+## SAVE:      Saves a copy of the problem if SAVE<>0.
+##            The file name can not be specified and defaults to "outpb.lp".
+##            The output file is CPLEX LP format.
+##
+## --- OUTPUTS ---
+##
+## XOPT:      The optimizer.
+##
+## FOPT:      The optimum.
+##
+## STATUS:    Status of the optimization.
+##     
+##              - Simplex Method -
+##              Value   Code     
+##              180     LPX_OPT     solution is optimal   
+##              181     LPX_FEAS    solution is feasible
+##              182     LPX_INFEAS  solution is infeasible
+##              183     LPX_NOFEAS  problem has no feasible solution
+##              184     LPX_UNBND   problem has no unbounded solution
+##              185     LPX_UNDEF   solution status is undefined
+##
+##              - Interior Point Method -
+##              Value   Code
+##              150     LPX_T_UNDEF the interior point method is undefined
+##              151     LPX_T_OPT   the interior point method is optimal
+##              *  Note that additional status codes may appear in
+##              the future versions of this routine *
+##    
+##              - Mixed Integer Method -
+##              Value   Code
+##              170     LPX_I_UNDEF  the status is undefined
+##              171     LPX_I_OPT    the solution is integer optimal
+##              172     LPX_I_FEAS   solution integer feasible but
+##                                   its optimality has not been proven
+##              173     LPX_I_NOFEAS no integer feasible solution
+##              
+## EXTRA:     A data structure containing the following fields:
+##             LAMBDA     Dual variables
+##             REDCOSTS   Reduced Costs
+##             TIME       Time (in seconds) used for solving LP/MIP problem in seconds.
+##             MEM        Memory (in bytes) used for solving LP/MIP problem.
+##
+##
+## In case of error the glpkmex returns one of the
+## following codes (these codes are in STATUS). For more informations on
+## the causes of these codes refer to the GLPK reference manual.
+##
+## Value  Code
+## 204    LPX_E_FAULT  unable to start the search
+## 205    LPX_E_OBJLL  objective function lower limit reached
+## 206    LPX_E_OBJUL  objective function upper limit reached
+## 207    LPX_E_ITLIM  iterations limit exhausted
+## 208    LPX_E_TMLIM  time limit exhausted
+## 209    LPX_E_NOFEAS no feasible solution
+## 210    LPX_E_INSTAB numerical instability
+## 211    LPX_E_SING   problems with basis matrix
+## 212    LPX_E_NOCONV no convergence (interior)
+## 213    LPX_E_NOPFS  no primal feas. sol. (LP presolver)
+## 214    LPX_E_NODFS  no dual feas. sol.   (LP presolver)
+##
+
+% --- CODE ---
+% If there is no input output the version and syntax
+if nargin==0
+	printf("glpk: An Octave interface to the GLPK library\n");
+	printf("Version: 1.0\n");
+	printf("\nSyntax: [xopt,fopt,status,extra]=glpk(sense,c,a,b,ctype,lb,ub,vartype,param,lpsolver,save)\n");
+	return;
+endif
+
+% If there are less than 4 arguments output an error message
+if nargin<4
+	error("At least 4 inputs required (sense,c,a,b)\n");
+	return;
+endif
+
+% At least 2 outputs required
+if nargout<2
+	error("2 outputs required\n");
+	return;
+endif
+
+% 1) Sense of optimization
+sense=varargin{1};
+
+% 2) Cost vector
+c=varargin{2};
+
+if (all(size(c)>1) | iscomplex(c) | ischar(c))
+	error("C must be a real vector\n");
+	return;
+endif
+nx=length(c);
+if size(c,1) ~= nx
+	c=c';
+endif
+
+% 3) Matrix constraint
+a=varargin{3};
+if isempty(a)
+	error("A cannot be an empty matrix\n");
+	return;
+endif
+[nc, nxa]=size(a);
+if (ischar(a) | iscomplex(a) | (nxa ~= nx))
+	error("A must be a real valued %d by %d matrix\n",nc,nx);
+	return;
+endif
+
+% 4) RHS
+b=varargin{4};
+if isempty(b)
+	error("B cannot be an empty vector\n");
+	return;
+endif
+if (ischar(b) | iscomplex(b) | (length(b) ~= nc))
+	error("B must be a real valued %d by 1 vector\n",nc);
+	return;
+endif
+
+% 5) Sense of each constraint
+ctype=[];
+if nargin>4 
+	ctype=varargin{5};
+	if isempty(ctype)
+		ctype=char('U'*ones(nc,1));
+	elseif (isnumeric(ctype) | all(size(ctype)>1) | (length(ctype) ~= nc))
+		error("CTYPE must be a char valued %d by 1 column vector\n",nc);
+		return;
+	else
+		for i=1:nc
+			if (ctype(i)!='F' & ctype(i)!='U' & ctype(i)!='S' & ctype(i)!='L' & ctype(i)!='D')
+				error("CTYPE must contain only F,U,S,L and D\n");
+				return;
+			endif
+		endfor
+	endif
+else
+	ctype=char('U'*ones(nc,1));
+end
+
+% 6) Vector with the lower bound of each variable
+lb=[];
+if nargin>5
+	lb=varargin{6};
+	if isempty(lb)
+		lb=-Inf*ones(nx,1);
+	elseif (ischar(lb) | iscomplex(lb) | all(size(lb)>1) | (length(lb)~=nx))
+		error("LB must be a real valued %d by 1 column vector\n",nx);
+		return;
+	endif
+else
+	lb=-Inf*ones(nx,1);
+end
+
+% 7) Vector with the upper bound of each variable
+ub=[];
+if nargin>6
+	ub=varargin{7};
+	if isempty(ub)
+		ub=Inf*ones(nx,1);
+	elseif (ischar(ub) | iscomplex(ub) | all(size(ub)>1) | (length(ub)~=nx))
+		error("UB must be a real valued %d by 1 column vector\n",nx);
+		return;
+	endif
+else
+	ub=Inf*ones(nx,1);
+end
+
+% 8) Vector with the type of variables 
+vartype=[];
+if nargin>7
+	vartype=varargin{8};
+	if isempty(vartype)
+		vartype=char('C'*ones(nx,1));
+	elseif (isnumeric(vartype) | all(size(vartype)>1) | (length(vartype)~=nx))
+		error("VARTYPE must be a char valued %d by 1 column vector\n",nx);
+		return;
+	else
+		for i=1:nx
+			if (vartype(i)!='C' & vartype(i)!='I')
+				error("VARTYPE must contain only C or I\n");
+				return;
+			endif
+		endfor
+	endif
+else
+	vartype=char('C'*ones(nx,1)); % As default we consider continuous vars
+endif
+
+% 9) Parameters vector
+param=[];
+if nargin>8
+	param=varargin{9};
+	if !isstruct(param)
+		error("PARAM must be a structure\n");
+		return;
+	endif
+else
+	param=struct;
+endif
+
+% 10) Select solver method: simplex or interior point
+lpsolver=[];
+if nargin>9
+	lpsolver=varargin{10};
+	if (!isnumeric(lpsolver) | all(size(lpsolver)>1))
+		error("LPSOLVER must be a real scalar value\n");
+		return;
+	endif
+else
+	lpsolver=1;
+endif
+
+% 11) Save the problem
+savepb=[];
+if nargin>10
+	savepb=varargin{11};
+	if (!isnumeric(savepb) | all(size(savepb)>1))
+		error("LPSOLVER must be a real scalar value\n");
+		return;
+	endif
+else
+	savepb=0;
+end
+
+if nargin>11
+	warning("Extra parameters ignored\n");
+endif
+
+try
+	[xopt, fmin, status, extra]=glpkoct(sense,c,a,b,ctype,lb,ub,vartype,param,lpsolver,savepb);
+catch
+	error("Problems with glpkoct");
+end_try_catch
+
+varargout{1}=xopt;
+varargout{2}=fmin;
+varargout{3}=status;
+varargout{4}=extra;
+
+	
+endfunction