Review or Modify Optimization Problems

Review Problem Using show or write
Change Default Solver or Options
Correct a Misspecified Problem
Duplicate Variable Name

Review Problem Using `show` or `write`

After you create an optimization problem, you can review its formulation by using show. For large problems, use write instead. For example,

prob = optimproblem;
x = optimvar('x',2,'LowerBound',0);
prob.Objective = x(1) - 2*x(2);
prob.Constraints.cons1 = x(1) + 2*x(2) <= 4;
prob.Constraints.cons2 = -x(1) + x(2) <= 1;

show(prob)

  OptimizationProblem : 

	Solve for:
       x

	minimize :
       x(1) - 2*x(2)


	subject to cons1:
       x(1) + 2*x(2) <= 4

	subject to cons2:
       -x(1) + x(2) <= 1

	variable bounds:
       0 <= x(1)
       0 <= x(2)

This review shows the basic elements of the problem, such as whether the problem is to minimize or maximize, and the variable bounds. The review shows the index names, if any, used in the variables. The review does not show whether the variables are integer valued.

Change Default Solver or Options

To attempt to improve a solution or speed of solution, examine and change the default solver or options.

To see the default solver and options, use optimoptions(prob). For example,

rng default
x = optimvar('x',3,'LowerBound',0);
expr = sum((rand(3,1).*x).^2);
prob = optimproblem('Objective',expr);
prob.Constraints.lincon = sum(sum(randn(size(x)).*x)) <= randn;
options = optimoptions(prob)

options = 

  lsqlin options:

   Options used by current Algorithm ('interior-point'):
   (Other available algorithms: 'trust-region-reflective')

   Set properties:
     No options set.

   Default properties:
              Algorithm: 'interior-point'
    ConstraintTolerance: 1.0000e-08
                Display: 'final'
           LinearSolver: 'auto'
          MaxIterations: 200
    OptimalityTolerance: 1.0000e-08
          StepTolerance: 1.0000e-12

   Show options not used by current Algorithm ('interior-point')

The default solver for this problem is lsqlin, and you can see the default options.

To change the solver, set the 'Solver' name-value pair in solve. To see the applicable options for a different solver, use optimoptions to pass the current options to the different solver. For example, continuing the problem,

options = optimoptions('quadprog',options)

options = 

  quadprog options:

   Options used by current Algorithm ('interior-point-convex'):
   (Other available algorithms: 'trust-region-reflective')

   Set properties:
    ConstraintTolerance: 1.0000e-08
          MaxIterations: 200
    OptimalityTolerance: 1.0000e-08
          StepTolerance: 1.0000e-12

   Default properties:
              Algorithm: 'interior-point-convex'
                Display: 'final'
           LinearSolver: 'auto'

   Show options not used by current Algorithm ('interior-point-convex')

To change the options, use optimoptions or dot notation to set options, and pass the options to solve in the 'Options' name-value pair. See Options in Common Use: Tuning and Troubleshooting. Continuing the example,

options.Display = 'iter';
sol = solve(prob,'Options',options,'Solver','quadprog');

 Iter            Fval  Primal Infeas    Dual Infeas  Complementarity  
    0    1.500359e+00   3.068423e-01   2.275437e+00     2.500000e-01  
    1    1.728717e-01   0.000000e+00   7.719860e-03     3.637874e-02  
    2    2.604108e-02   0.000000e+00   0.000000e+00     5.245260e-03  
    3    7.822161e-03   0.000000e+00   2.775558e-17     1.407915e-03  
    4    2.909218e-03   0.000000e+00   6.938894e-18     2.070784e-04  
    5    1.931264e-03   0.000000e+00   1.734723e-18     2.907724e-05  
    6    1.797508e-03   0.000000e+00   2.602085e-18     4.083167e-06  
    7    1.775398e-03   0.000000e+00   4.336809e-19     5.102453e-07  
    8    1.772971e-03   0.000000e+00   2.632684e-19     3.064243e-08  
    9    1.772848e-03   0.000000e+00   5.228973e-19     4.371356e-11  

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.

Correct a Misspecified Problem

To check that your problem is correct, review all its aspects. For example, create an optimization problem to solve a Sudoku problem by running this script.

x = optimvar('x',9,9,9,'LowerBound',0,'UpperBound',1);
cons1 = sum(x,1) == 1;
cons2 = sum(x,2) == 1;
cons3 = sum(x,3) == 1;
prob = optimproblem;
prob.Constraints.cons1 = cons1;
prob.Constraints.cons2 = cons2;
prob.Constraints.cons3 = cons3;
mul = ones(1,1,9);
mul = cumsum(mul,3);
prob.Objective = sum(sum(sum(x,1),2).*mul);
cons4 = optimconstr(3,3,9);

for u = 1:3
    for v = 1:3
        arr = x(3*(u-1)+1:3*(u-1)+3,3*(v-1)+1:3*(v-1)+3,:);
        cons4(u,v,:) = sum(sum(arr,1),2) <= ones(1,1,9);
    end
end
prob.Constraints.cons4 = cons4;

B = [1,2,2;
1,5,3;
1,8,4;
2,1,6;
2,9,3;
3,3,4;
3,7,5;
4,4,8;
4,6,6;
5,1,8;
5,5,1;
5,9,6;
6,4,7;
6,6,5;
7,3,7;
7,7,6;
8,1,4;
8,9,8;
9,2,3;
9,5,4;
9,8,2];

for u = 1:size(B,1)
    x.LowerBound(B(u,1),B(u,1),B(u,1)) = 1;
end

This script has some errors that you can find by examining the variables, objective, and constraints. First, examine the variable x.

x = 

  9×9×9 OptimizationVariable array with properties:

 Array-wide properties:
          Name: 'x'
          Type: 'continuous'
    IndexNames: {{}  {}  {}}

 Elementwise properties:
    LowerBound: [9×9×9 double]
    UpperBound: [9×9×9 double]

See variables with show.
See bounds with showbounds.

This display shows that the type of the variable is continuous. The variable should be integer valued. Change the type.

x.Type = 'integer'

x = 

  9×9×9 OptimizationVariable array with properties:

 Array-wide properties:
          Name: 'x'
          Type: 'integer'
    IndexNames: {{}  {}  {}}

 Elementwise properties:
    LowerBound: [9×9×9 double]
    UpperBound: [9×9×9 double]

See variables with show.
See bounds with showbounds.

Check the bounds. There should be 21 lower bounds with the value 1, one for each row of B. Because x is a large array, write the bounds to a file instead of displaying them at the command line.

writebounds(x,'xbounds.txt')

Search the file xbounds.txt for all instances of 1 <=. Only nine lower bounds having the value 1, in the variables x(1,1,1), x(2,2,2), …, x(9,9,9). To investigate this discrepancy, examine the code where you set the lower bounds:

for u = 1:size(B,1)
    x.LowerBound(B(u,1),B(u,1),B(u,1)) = 1;
end

The line inside the loop should say x.LowerBound(B(u,1),B(u,2),B(u,3)) = 1;. Reset all lower bounds to zero, then run the corrected code.

x.LowerBound = 0;
for u = 1:size(B,1)
    x.LowerBound(B(u,1),B(u,2),B(u,3)) = 1;
end
writebounds(x,'xbounds.txt')

xbounds.txt now has the correct number of lower bound entries that are 1.

Examine the objective function. The objective function expression is large, so write the expression to a file.

write(prob.Objective,'objectivedescription.txt')

    x(1, 1, 1) + x(2, 1, 1) + x(3, 1, 1) + x(4, 1, 1) + x(5, 1, 1) + x(6, 1, 1) + x(7, 1, 1) + x(8,
    1, 1) + x(9, 1, 1) + x(1, 2, 1) + x(2, 2, 1) + x(3, 2, 1) + x(4, 2, 1) + x(5, 2, 1) + x(6, 2,
    ...
    9*x(7, 8, 9) + 9*x(8, 8, 9) + 9*x(9, 8, 9) + 9*x(1, 9, 9) + 9*x(2, 9, 9) + 9*x(3, 9, 9) +
    9*x(4, 9, 9) + 9*x(5, 9, 9) + 9*x(6, 9, 9) + 9*x(7, 9, 9) + 9*x(8, 9, 9) + 9*x(9, 9, 9)

The objective function looks reasonable, because it is a sum of scalar expressions.

Write the constraints to files for examination.

write(prob.Constraints.cons1,'cons1.txt')
write(prob.Constraints.cons2,'cons2.txt')
write(prob.Constraints.cons3,'cons3.txt')
write(prob.Constraints.cons4,'cons4.txt')

Review cons4.txt and you see a mistake. All the constraints are inequalities rather than equalities. Correct the lines of code that create this constraint and put the corrected constraint in the problem.

cons4 = optimconstr(3,3,9);

for u = 1:3
    for v = 1:3
        arr = x(3*(u-1)+1:3*(u-1)+3,3*(v-1)+1:3*(v-1)+3,:);
        cons4(u,v,:) = sum(sum(arr,1),2) == ones(1,1,9);
    end
end
prob.Constraints.cons4 = cons4;

After these changes, you can successfully solve the problem.

sol = solve(prob);
x = round(sol.x);
y = ones(size(x));
for k = 2:9
    y(:,:,k) = k; % multiplier for each depth k
end
S = x.*y; % multiply each entry by its depth
S = sum(S,3); % S is 9-by-9 and holds the solved puzzle

drawSudoku(S)

Duplicate Variable Name

If you recreate a variable, but already have an expression that uses the old variable, then you can get errors when incorporating the expressions into a single problem. See Variables with Duplicate Names Disallowed.