How to solve a non-linear optimization problem with matrix and vector decision variables?

Hey folks. I'm somewhat new to optimizations and got stuck trying to solve an optimization problem where I have two decision variables: one is a matrix and the other is a vector. While my objective function is a rather simple linear cost function, and most of my constraints are also rather easy to handle, there is one somewhat lengthy non-linear constraint that seems to cause issues.

The problem looks as follows:

In the following, capital symbols are matrices, vectors are indicated by arrows, and everything else is a scalar. Dimensions of the vectors and matrix variables are as follows:

In my specific example, the values of L and N are 39 and 46, but in other examples they may go up to several hundreds.

Objective function: minimize total costs resulting from changes in b and p

Constraints:

1. Min and max allowed values for elements of DeltaB-matrix:
2. All elements outside of the diagonal of DeltaB-matrix must remain zero:
3. Min and max allowed values in elements of p-vector:
4. Changes in p-elements must sum up to zero:
5. And now the somewhat convoluted non-linear constraint which is rather difficult to explain:

Here's an MWE to showcase how I hoped to solve this using ga (I originally started with fmincon but that threw errors and recommended ga):

%% Input values
B = [ 10   -1  -2  -3  -4;...
       -1  19  -5  -6  -7;...
       -2  -5  24  -8  -9;...
       -3  -6  -8  27 -10;...
       -4  -7  -9 -10  30];
A = [  1  -1   0  0;...
       0   1  -1  0;...
       0   0   1 -1;...
       1   0   0 -1;...
       0   1   0 -1];
p = [-60; 160; -200; 100];
L = size(B,1); % = 5
N = size(A,2); % = 4
cost_per_b = 0.01; % same for all b_ij
delta_B_min = diag(-10*ones(L,1));
delta_B_max = diag(10*ones(L,1));
delta_p_min = -5*ones(N,1);
delta_p_max = 5*ones(N,1);
m = 120*ones(L,1);
cost_per_p = 1; % same for all elements of delta_p
mstart = -B * A * (A' * B * A)^-1 * p; % to give an idea for m-ranges

%% Optimization problem
% Decision variables
delta_B = optimvar('delta_B',L,L);
delta_p = optimvar('delta_p',N,1);
% Starting points
x0.delta_B = zeros(L);
x0.delta_p = zeros(N,1);
% Constraints
delta_B.LowerBound = delta_B_min;
delta_B.UpperBound = delta_B_max;
delta_p.LowerBound = delta_p_min;
delta_p.UpperBound = delta_p_max;
pbalance_cons = sum(delta_p) == 0;
nonlincons_pos = ...
    ((B + delta_B) * A * (A' * (B + delta_B) * A)^(-1)) * (p + delta_p)...
    <= m;
nonlincons_neg = ...
    ((B + delta_B) * A * (A' * (B + delta_B) * A)^(-1)) * (p + delta_p)...
    >= -m;
% Objective
objective = cost_per_b*sum((delta_B(:)')) + cost_per_p*sum(delta_p');

% Create optimization problem
prob = optimproblem('Objective',objective);
prob.Constraints.pbalance = pbalance_cons;
prob.Constraints.nonlincons_pos = nonlincons_pos;
prob.Constraints.nonlincons_neg = nonlincons_neg;

opts = optimoptions('ga');

show(prob);

% Solve optimization problem
[sol,fval,exitflag,output] = solve(prob,x0,opts);

While the output of show() looks fine (imo), the last line unfortunately produces an output error:

Error using optim.internal.problemdef.ProblemImpl/solveImpl

The value of 'GlobalSolver' is invalid. GlobalSolver must be a MultiStart or GlobalSearch object.

Error in optim.problemdef.OptimizationProblem/solve

Error in mwe_optimization (line 57)

[sol,fval,exitflag,output] = solve(prob,x0,opts);

I've tried to understand how to turn this into a MultiStart or GlobalSearch but I don't quite get what makes my optimization problem so difficult for just regular ga or other solvers. Why do I need MultiStart or GlobalSearch for this? And if it is really necessary, how do I set GlobalSolver accordingly?

Answer 1

[Torsten]

I turned deltaB into a diagonal matrix.

I changed m to 130 and it works. Are you sure a solution exists for m = 120 ?

%% Input values
B = [ 10   -1  -2  -3  -4;...
       -1  19  -5  -6  -7;...
       -2  -5  24  -8  -9;...
       -3  -6  -8  27 -10;...
       -4  -7  -9 -10  30];
A = [  1  -1   0  0;...
       0   1  -1  0;...
       0   0   1 -1;...
       1   0   0 -1;...
       0   1   0 -1];
p = [-60; 160; -200; 100];
L = size(B,1); % = 5
N = size(A,2); % = 4
cost_per_b = 0.01; % same for all b_ij
delta_B_min = -10*ones(L,1);
delta_B_max = 10*ones(L,1);
delta_p_min = -5*ones(N,1);
delta_p_max = 5*ones(N,1);
m = 130*ones(L,1);
cost_per_p = 1; % same for all elements of delta_p
mstart = -B * A * (A' * B * A)^-1 * p % to give an idea for m-ranges

mstart = 5×1

15 -130 45 5 -55

%% Optimization problem
% Decision variables
delta_B = optimvar('delta_B',L,1);
delta_p = optimvar('delta_p',N,1);
% Starting points
x0.delta_B = zeros(L,1);
x0.delta_p = zeros(N,1);
% Constraints
delta_B.LowerBound = delta_B_min;
delta_B.UpperBound = delta_B_max;
delta_p.LowerBound = delta_p_min;
delta_p.UpperBound = delta_p_max;
pbalance_cons = sum(delta_p) == 0;
nonlincons_pos = ...
    ((B + diag(delta_B)) * A * (A' * (B + diag(delta_B)) * A)^(-1)) * (p + delta_p)...
    <= m;
nonlincons_neg = ...
    ((B + diag(delta_B)) * A * (A' * (B + diag(delta_B)) * A)^(-1)) * (p + delta_p)...
    >= -m;
% Objective
objective = cost_per_b*sum(delta_B') + cost_per_p*sum(delta_p');

% Create optimization problem
prob = optimproblem('Objective',objective);
prob.Constraints.pbalance = pbalance_cons;
prob.Constraints.nonlincons_pos = nonlincons_pos;
prob.Constraints.nonlincons_neg = nonlincons_neg;

%opts = optimoptions('ga');

show(prob);

  OptimizationProblem : 

	Solve for:
       delta_B, delta_p

	minimize :
       0.01*delta_B(1) + 0.01*delta_B(2) + 0.01*delta_B(3) + 0.01*delta_B(4) + 0.01*delta_B(5) + delta_p(1) + delta_p(2) + delta_p(3) + delta_p(4)


	subject to pbalance:
       delta_p(1) + delta_p(2) + delta_p(3) + delta_p(4) == 0

	subject to nonlincons_pos:
       arg_LHS <= extraParams{7}

       where:

             arg1 = (extraParams{6} + delta_p);
             arg_LHS = ((((extraParams{4} + diag(delta_B)) * extraParams{5}) * ((extraParams{2} * (extraParams{1} + diag(delta_B))) * extraParams{3})^-1) * arg1);

         extraParams{1}:

         10    -1    -2    -3    -4    -1    19    -5    -6    -7    -2    -5    24    -8    -9    -3    -6    -8    27   -10

         :
         :

       extraParams{2}:

          1     0     0     1     0
         -1     1     0     0     1
          0    -1     1     0     0
          0     0    -1    -1    -1


       extraParams{3}:

          1    -1     0     0
          0     1    -1     0
          0     0     1    -1
          1     0     0    -1
          0     1     0    -1


       extraParams{4}:

         10    -1    -2    -3    -4    -1    19    -5    -6    -7    -2    -5    24    -8    -9    -3    -6    -8    27   -10

         :
         :

       extraParams{5}:

          1    -1     0     0
          0     1    -1     0
          0     0     1    -1
          1     0     0    -1
          0     1     0    -1


       extraParams{6}:

        -60
        160
       -200
        100


       extraParams{7}:

        130
        130
        130
        130
        130




	subject to nonlincons_neg:
       arg_LHS >= extraParams{7}

       where:

             arg1 = (extraParams{6} + delta_p);
             arg_LHS = ((((extraParams{4} + diag(delta_B)) * extraParams{5}) * ((extraParams{2} * (extraParams{1} + diag(delta_B))) * extraParams{3})^-1) * arg1);

         extraParams{1}:

         10    -1    -2    -3    -4    -1    19    -5    -6    -7    -2    -5    24    -8    -9    -3    -6    -8    27   -10

         :
         :

       extraParams{2}:

          1     0     0     1     0
         -1     1     0     0     1
          0    -1     1     0     0
          0     0    -1    -1    -1


       extraParams{3}:

          1    -1     0     0
          0     1    -1     0
          0     0     1    -1
          1     0     0    -1
          0     1     0    -1


       extraParams{4}:

         10    -1    -2    -3    -4    -1    19    -5    -6    -7    -2    -5    24    -8    -9    -3    -6    -8    27   -10

         :
         :

       extraParams{5}:

          1    -1     0     0
          0     1    -1     0
          0     0     1    -1
          1     0     0    -1
          0     1     0    -1


       extraParams{6}:

        -60
        160
       -200
        100


       extraParams{7}:

       -130
       -130
       -130
       -130
       -130




	variable bounds:
       -10 <= delta_B(1) <= 10
       -10 <= delta_B(2) <= 10
       -10 <= delta_B(3) <= 10
       -10 <= delta_B(4) <= 10
       -10 <= delta_B(5) <= 10

       -5 <= delta_p(1) <= 5
       -5 <= delta_p(2) <= 5
       -5 <= delta_p(3) <= 5
       -5 <= delta_p(4) <= 5

% Solve optimization problem
[sol,fval,exitflag,output] = solve(prob,x0);%,opts);

Solving problem using fmincon.

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

[Oliver Pohl]

Thank you, this actually fixed a problem I had encountered before when I tried to solve this with fmincon. I had tried the same thing but instead of using diag() I tried to manually turn an Lx1-vector of delta_b values into a diagonal matrix with

deltaB = eye(L) .* delta_b;

but fmincon appearantly doesn't like piecewise multiplications, so I switched to ga. With your solution I can go back to fmincon.

Unfortunately, it seems I didn't properly derive the MWE from my actual (a little bit more complex) code, because the MWE works fine now, but I still have issues with my actual code. I will work on it and update this.

[Oliver Pohl]

Ok, so I worked a bit more on this and got stuck on another problem: I got the non-linear constraints a bit wrong in my initial question, and the matrix inversion in them is actually not a "true" inversion but has to be a pseudo-inverse. Unfortunately, fmincon does not support the use of pinv(). Is there another solver that supports pseudo-inverted matrices?

I've also tried getting the inverted part containing the decision variable delta_b to the other side of the equation but I think I'm messing up the commutation laws because what I did on paper does not hold up when trying it with some examplary values in Matlab. If there is no solver capable of dealing with pseudo-inverted matrices, is there a way to remove the pseudo-inverse from my constraints?

To make things at least a little bit less convoluted, I added the dimensions to the matrices:

Here are the corrected non-linear constraints:

nonlincons_pos = ...
    ( (B + diag(delta_b)) * A * pinv(A' * (B + diag(delta_b)) * A) ) * (p + delta_p)...
    <= m;
nonlincons_neg = ...
    ( (B + diag(delta_b)) * A * pinv(A' * (B + diag(delta_b)) * A) ) * (p + delta_p)...
    >= -m;

This is how I tried removing the pinv() with delta_b in it:

Also, in my original question I got some variables' values wrong. Here are the corrected variables:

B = [ -10   0   0   0   0;...
        0 -11   0   0   0;...
        0   0 -12   0   0;...
        0   0   0 -13   0;...
        0   0   0   0 -14];

p = [1.32; -1.7; -2.0; 2.3];

% Give the elements 1,3,4,5 some room but restrict element 2
[~,m_max_idx] = max(abs(mstart));
m = mstart * 2;
m(m_max_idx) = 0.9*mstart(m_max_idx);

[Torsten]

%% Input values
B = [ -10   0   0   0   0;...
    0 -11   0   0   0;...
    0   0 -12   0   0;...
    0   0   0 -13   0;...
    0   0   0   0 -14];
p = [1.32; -1.7; -2.0; 2.3];
A = [  1  -1   0  0;...
    0   1  -1  0;...
    0   0   1 -1;...
    1   0   0 -1;...
    0   1   0 -1];
L = size(B,1); % = 5
N = size(A,2); % = 4
cost_per_b = 0.01; % same for all b_ij
delta_b_min = -10*ones(L,1);
delta_b_max = 10*ones(L,1);
delta_p_min = -5*ones(N,1);
delta_p_max = 5*ones(N,1);
m = 130*ones(L,1);
cost_per_p = 1; % same for all elements of delta_p
mstart = -B * A * (A' * B * A)^-1 * p % to give an idea for m-ranges

mstart = 5×1

-1.5850 -0.5000 2.4800 -1.8500 1.8100

% Give the elements 1,3,4,5 some room but restrict element 2
[~,m_max_idx] = max(abs(mstart));
m = mstart * 2;
m(m_max_idx) = 0.9*mstart(m_max_idx);
%% Optimization problem
% Decision variables
delta_b = optimvar('delta_b',L,1);
delta_p = optimvar('delta_p',N,1);
% Starting points
x0.delta_b = zeros(L,1);
x0.delta_p = zeros(N,1);
% Constraints
delta_b.LowerBound = delta_b_min;
delta_b.UpperBound = delta_b_max;
delta_p.LowerBound = delta_p_min;
delta_p.UpperBound = delta_p_max;
pbalance_cons = sum(delta_p) == 0;
inv(A'*(B+x0.delta_b)*A)

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.816655e-17.

ans = 4×4

1.0e+14 * -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643

rank(A'*(B+x0.delta_b)*A)

ans = 3

nonlincons = @(x,y)( (B + diag(x)) * A * pinv(A' * (B + diag(x)) * A) ) * (p + y);
nonlincons = fcn2optimexpr(nonlincons,delta_b,delta_p);
nonlincons_pos = nonlincons <= m;
nonlincons_neg = nonlincons >= -m;
% Objective
objective = cost_per_b*sum(delta_b') + cost_per_p*sum(delta_p');
% Create optimization problem
prob = optimproblem('Objective',objective);
prob.Constraints.pbalance = pbalance_cons;
prob.Constraints.nonlincons_pos = nonlincons_pos;
prob.Constraints.nonlincons_neg = nonlincons_neg;
%opts = optimoptions('ga');
show(prob);

  OptimizationProblem : 

	Solve for:
       delta_b, delta_p

	minimize :
       0.01*delta_b(1) + 0.01*delta_b(2) + 0.01*delta_b(3) + 0.01*delta_b(4) + 0.01*delta_b(5) + delta_p(1) + delta_p(2) + delta_p(3) + delta_p(4)


	subject to pbalance:
       delta_p(1) + delta_p(2) + delta_p(3) + delta_p(4) == 0

	subject to nonlincons_pos:
       arg_LHS <= extraParams{1}

       where:

             anonymousFunction1 = @(x,y)((B+diag(x))*A*pinv(A'*(B+diag(x))*A))*(p+y);
             arg_LHS = anonymousFunction1(delta_b, delta_p);

         extraParams{1}:

        -3.1700
        -1.0000
         2.2320
        -3.7000
         3.6200




	subject to nonlincons_neg:
       arg_LHS >= extraParams{1}

       where:

             anonymousFunction1 = @(x,y)((B+diag(x))*A*pinv(A'*(B+diag(x))*A))*(p+y);
             arg_LHS = anonymousFunction1(delta_b, delta_p);

         extraParams{1}:

         3.1700
         1.0000
        -2.2320
         3.7000
        -3.6200




	variable bounds:
       -10 <= delta_b(1) <= 10
       -10 <= delta_b(2) <= 10
       -10 <= delta_b(3) <= 10
       -10 <= delta_b(4) <= 10
       -10 <= delta_b(5) <= 10

       -5 <= delta_p(1) <= 5
       -5 <= delta_p(2) <= 5
       -5 <= delta_p(3) <= 5
       -5 <= delta_p(4) <= 5

% Solve optimization problem
[sol,fval,exitflag,output] = solve(prob,x0);%,opts);

Solving problem using fmincon.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.496562e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.995624e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.766987e-18.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.579935e-16.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.744714e-18.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.430329e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.636843e-19.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.051755e-17.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.641931e-19.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Warning: Matrix is singular to working precision.

Converged to an infeasible point.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance but constraints are not
satisfied to within the value of the constraint tolerance.


Consider enabling the interior point method feasibility mode.

[Oliver Pohl]

Thank you, this works! I can even add abs() into nonlcons so I don't have to explicitly put negative and positive restrictions separately:

m = abs(mstart) * 2; % should have actually done this from the start...
m(m_max_idx) = 0.9*abs(mstart(m_max_idx));
nonlincons = @(x,y)abs(((B + diag(x)) * A * pinv(A' * (B + diag(x)) * A) ) * (p + y));
nonlincons = fcn2optimexpr(nonlincons,delta_b,delta_p);
nonlincons_abs = nonlincons <= m;

With this, I can find a local minimum that satisfies the constraints for this MWE.

HOWEVER, in my real code, my matrices are larger than in the MWE, and here the solver stops prematurely:

Solving problem using fmincon.

Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,

options.MaxFunctionEvaluations = 3.000000e+03.

But when I try to increase MaxFunctionEvluations, I get an error:

opts = optimoptions('fmincon','MaxFunctionEvaluations',3e4,'MaxIterations',1500);

% Solve optimization problem
[sol,fval,exitflag,output] = solve(prob,x0,opts);

Error using optim.internal.problemdef.ProblemImpl/solveImpl

The value of 'GlobalSolver' is invalid. GlobalSolver must be a MultiStart or GlobalSearch object.

Error in optim.problemdef.OptimizationProblem/solve

Error in mwe_optimization (line 149)

[sol,fval,exitflag,output] = solve(prob,x0,opts);

This also happens in the MWE, and even when I only use opts = optimoptions('fmincon') without any other options. This is really confusing, because when I don't put any options, it uses fmincon anyways and it runs just fine. Why does this not work with options?

(Sorry for jumping from one problem to the next. I'm really new to this and trying my best. Thank you very much for your help!)

[Torsten]

Next time you ask a question please supply a complete executable code for which the problem appears so that we don't need to puzzle together the parts in order to test a possible solution.

%% Input values
B = [ -10   0   0   0   0;...
    0 -11   0   0   0;...
    0   0 -12   0   0;...
    0   0   0 -13   0;...
    0   0   0   0 -14];
p = [1.32; -1.7; -2.0; 2.3];
A = [  1  -1   0  0;...
    0   1  -1  0;...
    0   0   1 -1;...
    1   0   0 -1;...
    0   1   0 -1];
L = size(B,1); % = 5
N = size(A,2); % = 4
cost_per_b = 0.01; % same for all b_ij
delta_b_min = -10*ones(L,1);
delta_b_max = 10*ones(L,1);
delta_p_min = -5*ones(N,1);
delta_p_max = 5*ones(N,1);
cost_per_p = 1; % same for all elements of delta_p
mstart = -B * A * (A' * B * A)^-1 * p; % to give an idea for m-ranges
% Give the elements 1,3,4,5 some room but restrict element 2
[~,m_max_idx] = max(abs(mstart));
m = abs(mstart) * 2;
m(m_max_idx) = 0.9*abs(mstart(m_max_idx));
%% Optimization problem
% Decision variables
delta_b = optimvar('delta_b',L,1);
delta_p = optimvar('delta_p',N,1);
% Starting points
x0.delta_b = zeros(L,1);
x0.delta_p = zeros(N,1);
% Constraints
delta_b.LowerBound = delta_b_min;
delta_b.UpperBound = delta_b_max;
delta_p.LowerBound = delta_p_min;
delta_p.UpperBound = delta_p_max;
pbalance_cons = sum(delta_p) == 0;
inv(A'*(B+x0.delta_b)*A)

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.816655e-17.

ans = 4×4

1.0e+14 * -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643 -1.7643

rank(A'*(B+x0.delta_b)*A)

ans = 3

nonlincons = @(x,y)abs(((B + diag(x)) * A * pinv(A' * (B + diag(x)) * A) ) * (p + y));
nonlincons = fcn2optimexpr(nonlincons,delta_b,delta_p);
nonlincons_abs = nonlincons <= m;
% Objective
objective = cost_per_b*sum(delta_b') + cost_per_p*sum(delta_p');
% Create optimization problem
prob = optimproblem('Objective',objective);
prob.Constraints.pbalance = pbalance_cons;
prob.Constraints.nonlincons_pos = nonlincons_abs;
show(prob);

  OptimizationProblem : 

	Solve for:
       delta_b, delta_p

	minimize :
       0.01*delta_b(1) + 0.01*delta_b(2) + 0.01*delta_b(3) + 0.01*delta_b(4) + 0.01*delta_b(5) + delta_p(1) + delta_p(2) + delta_p(3) + delta_p(4)


	subject to pbalance:
       delta_p(1) + delta_p(2) + delta_p(3) + delta_p(4) == 0

	subject to nonlincons_pos:
       arg_LHS <= extraParams{1}

       where:

             anonymousFunction1 = @(x,y)abs(((B+diag(x))*A*pinv(A'*(B+diag(x))*A))*(p+y));
             arg_LHS = anonymousFunction1(delta_b, delta_p);

         extraParams{1}:

         3.1700
         1.0000
         2.2320
         3.7000
         3.6200




	variable bounds:
       -10 <= delta_b(1) <= 10
       -10 <= delta_b(2) <= 10
       -10 <= delta_b(3) <= 10
       -10 <= delta_b(4) <= 10
       -10 <= delta_b(5) <= 10

       -5 <= delta_p(1) <= 5
       -5 <= delta_p(2) <= 5
       -5 <= delta_p(3) <= 5
       -5 <= delta_p(4) <= 5

% Solve optimization problem
options = optimoptions(prob);
options.MaxFunctionEvaluations = 3e4;
options.MaxIterations = 1500;
[sol,fval,exitflag,output] = solve(prob,x0,'Options',options);

Solving problem using fmincon.

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

[Oliver Pohl]

Yes, you're right. I really did try to make my original example work though. Sorry for the back and forth.

Using options like that works now, thank you. Solver still exceeds the (increased) function evaluation limit, but I think my error lies elsewhere so I'll mark this as solved, since my original question has been thoroughly answered.