Optimization of a matrix with integers with nonlinear objective functions and nonlinear constraints

Hello all,

I have the following problem:

I want to optimize a matrix, which contains integer entries (e.g. A = [ 1 2 3 ; 2 1 1]), that the integers are arranged in the matrix. Each integer represents one area. I have linear and non-linear objective functions and boundary conditions.

My question now would be which algorithm would be useful. Intlinprog does not work in my opinion, because there are no non-linear objective functions or boundary conditions possible and fmincon does not work because you can not calculate any integer as x0. Do you have an idea which algorithm I can use in Matlab?

Thank you very much.

[Dyuman Joshi]

"fmincon does not work because you can not calculate any integer as x0."

I am not sure what you mean by this. Could you elaborate more?

[b999]

As a result I would like to get a matrix at the end, which contains as entries Integers. So for example the matrix A = [ 1 2 3 ; 2 1 1] shall be optimized to A = [ 2 2 1 ; 2 1 2]. And according to my research it is not possible with fmincon to get integers as final output. I understood it in such a way that fmincon would then output for example A = [ 2.1213 2.234526 1.151266 ; 2.165256 1.1416783 2.1515414]. However, I need A = [ 2 2 1 ; 2 1 2].

Answer 1

[John D'Errico]

You are correct, in that fmincon does not provide integer constraints as an option, and that a tool like intlinprog does not allow nonlinear constraints.

This class of problem will typically fall into the domain of tools from the Global Optimization toolbox, mainly GA. (I cannot recall if any other of those tools allow integer constraints. I checked once, and there were no other options at that time, except for GA.) GA can accept integer constraints on the variables, but also nonlinear constraints on them.

[Dyuman Joshi]

I found a few things by looking here and there (might not have been a thorough search)

surrogateopt

https://in.mathworks.com/help/gads/nonlinear-problem-integer-nonlinear-constraints.html#SolveNonlinearProblemWithIntegerAndNonlinearConstraintsExample-4

https://in.mathworks.com/matlabcentral/answers/845715-integer-based-optimization-solver-with-non-linear-constraints

Might be helpful.

Edit - https://in.mathworks.com/matlabcentral/answers/52443-force-fmincon-to-output-interger-values#answer_63928

[b999]

To John D'Errico

Thank you for your answer. By GA, do you mean the genetic algorithm?

I have two questions, if the genetic algorithm is meant:

1) I was wondering about this statement in GA: "When there are integer constraints, ga does not accept nonlinear equality constraints, only nonlinear inequality constraints". Is it then not possible to insert nonlinear equality constraints? I would probably need them.

2) Can a matrix be used as input or is it possible to get an optimized matrix as output? (as explained in the question).

Thank you very much.

[b999]

@Dyuman Joshi

surrogateopt sounds interesting. Never heard of it before. I am also relatively new to Matlab. Based on the description at "Help", which is

x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,intcon,options) for surrogateopt

reads, I have a hard time seeing if a matrix is possible as input/output. Often in this case we talk about vectors

[John D'Errico]

@a9v9

Yes, by GA, I meant the function named GA, which is a genetic algorithm.

You did not say which kind of nonlinear constraints. What you should understand is that nonlinear EQUALITY constraints where integers are concerned, are rather difficult to deal with.

For example, suppose your problem was to solve for some variables x and y (for a fixed value of n), such that

x^2 - n*y^2 == 1

This is a variation of Pell equation.

https://en.wikipedia.org/wiki/Pell%27s_equation

And depending on the value of n, just finding any feasible points can be not at all trivial to identify. It makes the problem MUCH more difficult. And so, if you have integer variables, then you cannot also have nonlinear EQUALITY constraints. You will find that to be a common issue.

And, yes, you might wonder why that is a problem at all. But my point is, for such a solver to work, it first needs to find valid points to evaluate the function at. And, if the problem was as I described, where the only valid feasible points are the solutions of a Pell equation, that makes it a very difficult problem. Sorry.

Can a matrix be used as input? Well, GA does not take starting values, as it is not that class of solver that needs them. But even if you could only get a solution in the form of a vector, then you could still trivially just reshape the vector into a matrix of your desired size and shape.

help ga

 GA    Constrained optimization using genetic algorithm.
    GA attempts to solve problems of the following forms:
        min F(X)  subject to:  A*X  <= B, Aeq*X  = Beq (linear constraints)
         X                     C(X) <= 0, Ceq(X) = 0 (nonlinear constraints)
                               LB <= X <= UB 
                               X(i) integer, where i is in the index
                               vector INTCON (integer constraints)
 
    Note: If INTCON is not empty, then no equality constraints are allowed.
    That is:-
    * Aeq and Beq must be empty
    * Ceq returned from NONLCON must be empty
 
    X = GA(FITNESSFCN,NVARS) finds a local unconstrained minimum X to the
    FITNESSFCN using GA. NVARS is the dimension (number of design
    variables) of the FITNESSFCN. FITNESSFCN accepts a vector X of size
    1-by-NVARS, and returns a scalar evaluated at X.
 
    X = GA(FITNESSFCN,NVARS,A,b) finds a local minimum X to the function
    FITNESSFCN, subject to the linear inequalities A*X <= B. Linear
    constraints are not satisfied when the PopulationType option is set to
    'bitString' or 'custom'. See the documentation for details.
 
    X = GA(FITNESSFCN,NVARS,A,b,Aeq,beq) finds a local minimum X to the
    function FITNESSFCN, subject to the linear equalities Aeq*X = beq as
    well as A*X <= B. (Set A=[] and B=[] if no inequalities exist.) Linear
    constraints are not satisfied when the PopulationType option is set to
    'bitString' or 'custom'. See the documentation for details.
 
    X = GA(FITNESSFCN,NVARS,A,b,Aeq,beq,lb,ub) defines a set of lower and
    upper bounds on the design variables, X, so that a solution is found in
    the range lb <= X <= ub. Use empty matrices for lb and ub if no bounds
    exist. Set lb(i) = -Inf if X(i) is unbounded below;  set ub(i) = Inf if
    X(i) is unbounded above. Linear constraints are not satisfied when the
    PopulationType option is set to 'bitString' or 'custom'. See the 
    documentation for details.
 
    X = GA(FITNESSFCN,NVARS,A,b,Aeq,beq,lb,ub,NONLCON) subjects the
    minimization to the constraints defined in NONLCON. The function
    NONLCON accepts X and returns the vectors C and Ceq, representing the
    nonlinear inequalities and equalities respectively. GA minimizes
    FITNESSFCN such that C(X)<=0 and Ceq(X)=0. (Set lb=[] and/or ub=[] if
    no bounds exist.) Nonlinear constraints are not satisfied when the
    PopulationType option is set to 'bitString' or 'custom'. See the 
    documentation for details.
 
    X = GA(FITNESSFCN,NVARS,A,b,Aeq,beq,lb,ub,NONLCON,options) minimizes
    with the default optimization parameters replaced by values in OPTIONS.
    OPTIONS can be created with the OPTIMOPTIONS function. See OPTIMOPTIONS
    for details. For a list of options accepted by GA refer to the
    documentation.
 
    X = GA(FITNESSFCN,NVARS,A,b,[],[],lb,ub,NONLCON,INTCON) requires that
    the variables listed in INTCON take integer values. Note that GA does
    not solve problems with integer and equality constraints. Pass empty
    matrices for the Aeq and beq inputs if INTCON is not empty.
 
    X = GA(FITNESSFCN,NVARS,A,b,[],[],lb,ub,NONLCON,INTCON,options)
    minimizes with integer constraints and the default optimization
    parameters replaced by values in OPTIONS. OPTIONS can be created with
    the OPTIMOPTIONS function. See OPTIMOPTIONS for details.
 
    X = GA(PROBLEM) finds the minimum for PROBLEM. PROBLEM is a structure
    that has the following fields:
        fitnessfcn: <Fitness function>
             nvars: <Number of design variables>
             Aineq: <A matrix for inequality constraints>
             bineq: <b vector for inequality constraints>
               Aeq: <Aeq matrix for equality constraints>
               beq: <beq vector for equality constraints>
                lb: <Lower bound on X>
                ub: <Upper bound on X>
           nonlcon: <Nonlinear constraint function>
            intcon: <Index vector for integer variables>
           options: <Options created with optimoptions('ga',...)>
          rngstate: <State of the random number generator>
 
    [X,FVAL] = GA(FITNESSFCN, ...) returns FVAL, the value of the fitness
    function FITNESSFCN at the solution X.
 
    [X,FVAL,EXITFLAG] = GA(FITNESSFCN, ...) returns EXITFLAG which
    describes the exit condition of GA. Possible values of EXITFLAG and the
    corresponding exit conditions are
 
      1 Average change in value of the fitness function over
         options.MaxStallGenerations generations less than 
         options.FunctionTolerance and constraint violation less than 
         options.ConstraintTolerance.
      3 The value of the fitness function did not change in
         options.MaxStallGenerations generations and constraint violation 
         less than options.ConstraintTolerance.
      4 Magnitude of step smaller than machine precision and constraint
         violation less than options.ConstraintTolerance. This exit 
         condition applies only to nonlinear constraints.
      5 Fitness limit reached and constraint violation less than
         options.ConstraintTolerance. 
      0 Maximum number of generations exceeded.
     -1 Optimization terminated by the output or plot function.
     -2 No feasible point found.
     -4 Stall time limit exceeded.
     -5 Time limit exceeded.
 
    [X,FVAL,EXITFLAG,OUTPUT] = GA(FITNESSFCN, ...) returns a
    structure OUTPUT with the following information:
              rngstate: <State of the random number generator before GA started>
           generations: <Total generations, excluding HybridFcn iterations>
             funccount: <Total function evaluations>
         maxconstraint: <Maximum constraint violation>, if any
               message: <GA termination message>
 
    [X,FVAL,EXITFLAG,OUTPUT,POPULATION] = GA(FITNESSFCN, ...) returns the
    final POPULATION at termination.
 
    [X,FVAL,EXITFLAG,OUTPUT,POPULATION,SCORES] = GA(FITNESSFCN, ...) returns
    the SCORES of the final POPULATION.
 
 
    Example:
      Unconstrained minimization of Rastrigins function:
        function scores = myRastriginsFcn(pop)
            scores = 10.0 * size(pop,2) + sum(pop.^2 - 10.0*cos(2*pi .* pop),2);
 
      numberOfVariables = 2
      x = ga(@myRastriginsFcn,numberOfVariables)
 
      Display plotting functions while GA minimizes
       options = optimoptions('ga','PlotFcn',...
         {@gaplotbestf,@gaplotbestindiv,@gaplotexpectation,@gaplotstopping});
       [x,fval,exitflag,output] = ga(fitfcn,2,[],[],[],[],[],[],[],options)
 
    An example with inequality constraints and lower bounds
     A = [1 1; -1 2; 2 1];  b = [2; 2; 3];  lb = zeros(2,1);
     fitfcn = @(x)0.5*x(1)^2 + x(2)^2 -x(1)*x(2) -2*x(1) - 6.0*x(2);
 
     % Use mutation function which can handle constraints
     options = optimoptions('ga','MutationFcn',@mutationadaptfeasible);
     [x,fval,exitflag] = ga(fitfcn,2,A,b,[],[],lb,[],[],options);
 
    If FITNESSFCN or NONLCON are parameterized, you can use anonymous
    functions to capture the problem-dependent parameters. Suppose you want
    to minimize the fitness given in the function myfit, subject to the
    nonlinear constraint myconstr, where these two functions are
    parameterized by their second argument a1 and a2, respectively. Here
    myfit and myconstr are MATLAB file functions such as
 
         function f = myfit(x,a1)
         f = exp(x(1))*(4*x(1)^2 + 2*x(2)^2 + 4*x(1)*x(2) + 2*x(2) + a1);
 
    and
 
         function [c,ceq] = myconstr(x,a2)
         c = [1.5 + x(1)*x(2) - x(1) - x(2);
               -x(1)*x(2) - a2];
         % No nonlinear equality constraints:
          ceq = [];
 
    To optimize for specific values of a1 and a2, first assign the values
    to these two parameters. Then create two one-argument anonymous
    functions that capture the values of a1 and a2, and call myfit and
    myconstr with two arguments. Finally, pass these anonymous functions to
    GA:
 
      a1 = 1; a2 = 10; % define parameters first
      % Mutation function for constrained minimization
      options = optimoptions('ga','MutationFcn',@mutationadaptfeasible);
      x = ga(@(x)myfit(x,a1),2,[],[],[],[],[],[],@(x)myconstr(x,a2),options)
 
    Example: Solving a mixed-integer optimization problem
    An example of optimizing a function where a subset of the variables are
    required to be integers:
    
    % Define the objective and call GA. Here variables x(2) and x(3) will
    % be integer. 
    fun = @(x) (x(1) - 0.2)^2 + (x(2) - 1.7)^2 + (x(3) -5.1)^2; 
    x = ga(fun,3,[],[],[],[],[],[],[],[2 3])
           
    See also OPTIMOPTIONS, FITNESSFUNCTION, GAOUTPUTFCNTEMPLATE, PATTERNSEARCH, @.

    Documentation for ga
       doc ga

So the fitness function is presumed to be a VECTOR, of length 1xnvars. But nothing stops you from reshaping that vector inside your objective function, and after it is returned from GA.

And, yes, surrogateopt does also allow integer constraints. However, it is designed a little differently, in that it has a target of problems with time-consuming functions. So if your function is not costly to evaluate, you may fund that GA does a little better at truly finding the global optimizer. I don't have any experience with surrogateopt, so I cannot know for sure. You might be well off to test your problem on both optimizers. I'd start with GA, unless your function is a costly one to evaluate.

help surrogateopt

 SURROGATEOPT The surrogateopt function is a global solver for
    time-consuming objective functions. surrogateopt attempts to solve
    problems of the form
 
    minimize f(x)
    subject to: lb <= x <= ub (bound constraints)
              A*X  <= B, Aeq*X  = Beq (linear constraints)
              c(x) <= 0 (nonlinear inequalities)
              x(i) must be integer for i in intcon
 
    The solver searches for the global minimum of a real-valued objective
    function in multiple dimensions, subject to bounds, optional integer
    constraints, and optional nonlinear inequality constraints. surrogateopt
    is best suited to objective functions that take a long time to evaluate.
    The objective function can be nonsmooth. The solver requires finite
    bounds on all variables. The solver can optionally maintain a checkpoint
    file to enable recovery from crashes or partial execution, or
    optimization continuation after meeting a stopping condition.
 
    x = SURROGATEOPT(fun,lb,ub) searches for a global minimum of fun(x) in
    the region lb <= x <= ub.
 
    fun accepts a single row vector argument x. fun can return either
     * a real scalar fval = fun(x), or
     * a structure. If the structure contains a field Fval, then
       surrogateopt attempts to minimize fun(x).Fval. If the structure
       contains a field Ineq, then surrogateopt attempts to make all of the
       components of that field be nonpositive: fun(x).Ineq <= 0 for all
       entries.
 
    For information on converting nonlinear constraints between the
    surrogateopt structure syntax and other solvers, see the packfcn
    function reference page.
 
    x = SURROGATEOPT(fun,lb,ub,intcon) requires that the variables listed
    in intcon take integer values. pass intcon = [] if there are no integer
    variables.
 
    x = SURROGATEOPT(fun,lb,ub,intcon,A,B) finds a local minimum x
    subject to the linear inequalities A*X <= B.
 
    x = SURROGATEOPT(fun,lb,ub,intcon,A,B,Aeq,Beq) finds a local minimum x
    subject to the linear equalities Aeq*X = Beq as well as A*X <= B. (Set
    A=[] and B=[] if no inequalities exist.)
 
    x = SURROGATEOPT(fun,lb,ub,intcon,A,B,Aeq,Beq,options) replaces
    the default optimization parameters by values in options. Create
    options using optimoptions.
 
    x = SURROGATEOPT(PROBLEM) searches for a minimum for PROBLEM, a
    structure with these fields: 
        objective: The objective (and constraint) function
               lb: Lower bounds for x
               ub: Upper bounds for x
           intcon: Indices of variables that must be integer
           Aineq: <A matrix for inequality constraints>
           bineq: <B vector for inequality constraints>
             Aeq: <A matrix for equality constraints>
             beq: <B vector for equality constraints>
        rngstate: Optional field to reset the state of RNG
         options: Options created with OPTIMOPTIONS
          solver: 'surrogateopt'
 
    x = SURROGATEOPT(checkpointFile) continues running the optimization
    from the state in a saved checkpoint file. checkpointFile is a path to
    a checkpoint file, specified as a string or character vector. If you
    specify a file name without a path, SURROGATEOPT uses a checkpoint file
    in the current folder.
 
    SURROGATEOPT creates a checkpoint file when it has a valid
    CheckpointFile option.
 
    The data in a checkpoint file is in .mat format. To avoid errors or
    other unexpected results, do not modify the data before calling
    surrogateopt.
 
    x = SURROGATEOPT(checkpointFile,opts) continues running the
    optimization from the state in a saved checkpoint file, and replaces
    options in checkpointFile with those in opts. surrogateopt accepts
    changes in the following options:
    CheckpointFile
    Display
    MaxFunctionEvaluations
    MaxTime
    MinSurrogatePoints
    ObjectiveLimit
    OutputFcn
    PlotFcn
    UseParallel
    UseVectorized
    BatchUpdateInterval
 
    [x,fval] = SURROGATEOPT(__) also returns the best (smallest) value of
    the objective function found by the solver, using any of the input
    argument combinations in the previous syntaxes.
 
    [x,fval,exitflag] = SURROGATEOPT(__) also returns exitflag, an
    integer describing the reason the solver stopped. Possible values of
    exitflag and the corresponding exit conditions are listed below.
 
      1  Best function value reached options.ObjectiveLimit
 
      3  Feasible point but solver stopped (stalled) because no new points
      were found.
      
      10  Unique solution exist for the problem due to one of the following
      reasons: 
         All lower bounds are same as corresponding upper bounds.
         Unique solution due to linear equalities.
 
      0  Number of function evaluations exceeded options.MaxFunctionEvaluations
         or elapsed time exceeded options.MaxTime.
 
     -1  Stopped by output/plot function.
 
     -2  No feasible point found due to one or more of the following reasons:
         One or more lower bound lb(i) exceeds a corresponding upper bound ub(i).
         One or more ceil(lb(i)) exceeds a corresponding floor(ub(i)) for i in INTCON.
         Linear constraints are infeasible.
         No solution found that satisfies nonlinear constraints.
 
    [x,fval,exitflag,output] = SURROGATEOPT(__) also returns output, a
    structure with these fields:
           rngstate: State of the random number generator before the solver started.
          funccount: Total number of function evaluations.
            message: Reason why SURROGATEOPT stopped.
        elapsedtime: Time spent running the solver in seconds, as measured by tic/toc.
    constrviolation: Maximum constraint violation due to nonlinear
                     inequalities at x.
               ineq: Inequality constraints at the x.
 
    [x,fval,exitflag,output,trials] = SURROGATEOPT(__) also returns a
    structure containing all of the evaluated points (Npts) and the
    function function values at those points. trials have these fields:
 
    trials.X: Matrix of size Npts-by-Nvar.
              Each row of trials.X represents one point that SURROGATEOPT
              evaluated.
    trials.Fval: A vector, where each entry is the objective function
                 value of the corresponding row of trials.X
    trials.Ineq: A matrix in which rows are nonlinear inequalties of the
                 corresponding row of trials.X
 
    Examples:
    Find the minimum of a function subject to bounds.
 
    Search for a minimum of the six-hump camel back function in the region
    -2.1 <= x(i) <= 2.1. This function has two global minima with the
    objective function value -1.0316284... and four local minima with
    higher objective function values.
 
    fun = @(x)(4*x(:,1).^2 - 2.1*x(:,1).^4 + x(:,1).^6/3 ...
               + x(:,1).*x(:,2) - 4*x(:,2).^2 + 4*x(:,2).^4);
 
    lb = [-2.1,-2.1];
    ub = -lb;
 
    x = surrogateopt(fun,lb,ub)
 
    Find the minimum of a function subject to nonlinear constraints.
 
    Find the minimum of Rosenbrock's function 100(x(2)-x(1)2)2+(1-x(1))2
    subject to the nonlinear constraint that the solution lies in a disk of
    radius 1/3 around the point [1/3,1/3]: (x(1)-1/3)2+(x(2)-1/3)2 <= (1/3)2.
 
    To do so, write a function fun(x) that returns the value of
    Rosenbrock's function in a structure field Fval, and returns the
    nonlinear constraint value in the form c(x) <= 0 in the structure field
    Ineq.
 
    function f = fun(x)
    f.Fval = 100*(x(2) - x(1)^2)^2 + (1 - x(1))^2;
    f.Ineq = (x(1)-1/3)^2 + (x(2)-1/3)^2 - (1/3)^2;
 
    Call surrogateopt using lower bounds of 0 on each component and upper
    bounds of 2/3.
    lb = [0,0];
    ub = [2/3,2/3];
 
    [x,fval,exitflag] = surrogateopt(@fun,lb,ub)
 
    If objective and constraints are in separate function use this
    convenienece function to return output as a struct.
 
    objAndConstr = packfcn(@objective,@inequalities)
 
    X = surrogateopt(objAndConstr, ...) % objAndConstr is a function handle
 
    Find the minimum of a function subject to integer constraints
 
    Find the minimum of the piecewiseFcn function for a two-dimensional
    variable x whose first component is restricted to integer values, and
    all components are between -5 and 5. 
  
          function f = piecewiseFcn(x)
          f = zeros(1,size(x,1));
          for i = 1:size(x,1)
              if  x(i,1) < -5
                  f(i) = (x(i,1)+5)^2 + abs(x(i,2));
              elseif x(i,1) < -3
                  f(i) = -2*sin(x(i,1)) + abs(x(i,2));
              elseif x(i,1) < 0
                  f(i) = 0.5*x(i,1) + 2 + abs(x(i,2));
              elseif x(i,1) >= 0
                  f(i) = .3*sqrt(x(i,1)) + 5/2 +abs(x(i,2));
              end
          end 
 
 
    intcon = 1;
    rng default % For reproducibility 
    fun = @piecewiseFcn;
    lb = [-5,-5];
    ub = [5,5];
    x = surrogateopt(fun,lb,ub,intcon)
  
    Restart from checkpoint file
  
    To enable restarting surrogate optimization due to a crash or any other
    reason, set a checkpoint file name.
 
    Create an optimization problem and set a small number of function
    evaluations.
    rng default % For reproducibility
  
    fun = @(x)(4*x(:,1).^2 - 2.1*x(:,1).^4 + x(:,1).^6/3 ...
            + x(:,1).*x(:,2) - 4*x(:,2).^2 + 4*x(:,2).^4);
    lb = [-2.1,-2.1];
    ub = -lb;
    opts = optimoptions('surrogateopt','CheckpointFile','checkfile.mat');
    opts.MaxFunctionEvaluations = 30;
  
    [x,fval,exitflag,output] = surrogateopt(fun,lb,ub,opts)
  
    Set options to use 100 function evaluations (which means 70 more than
    already done) and restart the optimization.
    opts.MaxFunctionEvaluations = 100;
 
    [x2,fval2,exitflag2,output2] = surrogateopt('checkfile.mat',opts)
 
    See also OPTIMOPTIONS, PATTERNSEARCH, FMINCON

    Documentation for surrogateopt
       doc surrogateopt

[Walter Roberson]

ga accepts an options structure that supports an 'InitialPopulationMatrix' -- so it can take starting values (just not conveniently)

The nonlinear equality A == B can be coded as the pair of nonlinear inequalities [A - B, B - A] which corresponds to A <= B & B <= A which is only true for A == B . You can do this -- but it will probably still have a lot of trouble finding matching points.