Finding solution for array that satisfies given constraints and gives desired output

Hello everyone,

the problem I am dealing with is that I have an array with versines measured:

EV = [0 2 -2 6 2 8 4 8 6 4 8 10 10 6 8 10 4 18 0 6 4 8 9 4 -4 6 16 4 -2 -2 -4 -6 -14 1 2 -2 -12 -8 -10 -14 -14 -18];

and I need to optimize this array to produce a smoother result with less spikes which would be a second array. However, this second array needs to have certain rules and just randomly making it smoother doesn't work because the following equations govern other parameters that are calculated from this:

d(1,1) = 0;
ED(1,1) = 0;
M(1,1) = 0;
S(1,1) = 0;
d(1,i) = DV(1,i)-EV(1,i);
ED(1,i) = ED(1,i-1) + d(1,i-1);
M(1,i) = ED(1,i) + M(1,i-1);
S(1,i) = -2*(ED(1,i-1)+DV(1,i-1)-EV(1,i-1)+M(1,i-1));

DV is the second (optimized/smoother) array that I need to create. The results of individual parameters are dependent on previous values. So for example:

ED(1,2) = ED(1,1) + d(1,1) 
M(1,2) = ED(1,2) + M(1,1)
S(1,2) = -2*(ED(1,1)+DV(1,1)-EV(1,1)+M(1,1))

And I need to create this second array such that the last two values of S wil be 0 so:

S(1,42) = 0 && S(1,41) = 0;

My approach:

I first created DV from a 6-point moving average of EV which provided a pretty good base and S at the end was -82. Then I wanted to compute all permutations close to this array by either adding 1, subtracting 1 on not changing the values. I believed this would tweak the result and I could then filter for the results which give two 0s at the end. However, with 42 positions and 3 possible values, this is simply not feasible and so I chose 12 positions in the middle of the array and 3 values. This created about a million permutations and some of them which were pretty close to what I need, my S at the end was 0 and my second to last S was 6. However, I cannot compute more permutations due to memory and time constraint. Is there any other approach like opitmization I can use to solve this?

Summary:

Create a second array that gives result of S at the last and second to last position 0.

Thank you for your help

Answer 1

[Alan Weiss]

I don't know if this is what you want, but the problem-based approach in Optimization Toolbox™ handles your equations easily.

EV = [0 2 -2 6 2 8 4 8 6 4 8 10 10 6 8 10 4 18 0 6 4 8 9 4 -4 6 16 4 -2 -2 -4 -6 -14 1 2 -2 -12 -8 -10 -14 -14 -18]';
% EV is a column to match the rest of the variables
% Create optimization variables
N = length(EV);
d = optimvar('d',N);
ED = optimvar('ED',N);
M = optimvar('M',N);
S = optimvar('S',N);
DV = optimvar('DV',N);
% Create constraints
cons1 = d == DV - EV;
consd0 = d(1) == 0;
consED0 = ED(1) == 0;
consM0 = M(1) == 0;
consS0 = S(1) == 0;
idxnot1 = 2:N;
consED = ED(idxnot1) == ED(idxnot1 - 1) + d(idxnot1 - 1);
consMnot1 = M(idxnot1) == ED(idxnot1) + M(idxnot1 - 1);
consSnot1 = S(idxnot1) == -2*(ED(idxnot1 - 1) + DV(idxnot1 - 1) - EV(idxnot1 - 1) + M(idxnot1 - 1));
% Create problem and put constraints in problem
prob = eqnproblem;
prob.Equations.cons1 = cons1;
prob.Equations.consd0 = consd0;
prob.Equations.consED0 = consED0;
prob.Equations.M0 = consM0;
prob.Equations.S0 = consS0;
prob.Equations.consED = consED;
prob.Equations.consMnot1 = consMnot1;
prob.Equations.consSnot1 = consSnot1;
prob.Equations.consSend = S(end-1:end) == [0;0];
% Solve problem
[sol,fval,eflag,output] = solve(prob);

Linear equation problem is underdetermined. Equations will be solved in a least squares sense.

Solving problem using lsqlin.

% View DV variables
disp(sol.DV)

     0
     2
    -2
     6
     2
     8
     4
     8
     6
     4
     8
    10
    10
     6
     8
    10
     4
    18
     0
     6
     4
     8
     9
     4
    -4
     6
    16
     4
    -2
    -2
    -4
    -6
   -14
     1
     2
    -2
   -12
   -22
     4
     0
   -28
   -18

Alan Weiss

MATLAB mathematical toolbox documentation

[Kevin Syc]

Hi Alan,

thank you so much for your solution, this is exactly the approach I am looking for, however the result it gives is not very smooth but rather very similar to the original line:

I would need the result look something like this:

Which still gives 0 for the last two digits of the slue and the orange line is much smoother. Is there a way to do it? Maybe with an initial guess of the 6 point average or something like that?

[Alan Weiss]

I don't understand what you are plotting. How do your plots relate to the problem variables?

In any case, you can feel free to modify the equations to get a different answer. You see that the original equations lead to an underdetermined system, meaning there is room for more constraints.

Alan Weiss

MATLAB mathematical toolbox documentation

[Kevin Syc]

Hi Alan,

sorry for the confusion, The blue line of this graph is the original measured versine (EV) and the orange one is the solution that your code provided. As you can see they are mostly the same as they are overlaying each other for the most part.

Instead, I would need the result to look something like this: Where the blue line is again the measured versine (EV) and the orange line is the result that your code would output. I would need it to be way smoother to have nice transition of the curve and remove any spikes, but still give me two 0s at the end of the Slue array.

In any way, thank you very much for your help.

[Alan Weiss]

Well, you can try entering the equations as constraints and solving a least-squares problem. One component of the least squares is the difference between DV and EV. The other is the difference between DV(j) and DV(j - 1). See if you like the result.

EV = [0 2 -2 6 2 8 4 8 6 4 8 10 10 6 8 10 4 18 0 6 4 8 9 4 -4 6 16 4 -2 -2 -4 -6 -14 1 2 -2 -12 -8 -10 -14 -14 -18]';
N = length(EV);
d = optimvar('d',N);
ED = optimvar('ED',N);
M = optimvar('M',N);
S = optimvar('S',N);
DV = optimvar('DV',N);

cons1 = d == DV - EV;
consd0 = d(1) == 0;
consED0 = ED(1) == 0;
consM0 = M(1) == 0;
consS0 = S(1) == 0;
idxnot1 = 2:N;
consED = ED(idxnot1) == ED(idxnot1 - 1) + d(idxnot1 - 1);
consMnot1 = M(idxnot1) == ED(idxnot1) + M(idxnot1 - 1);
consSnot1 = S(idxnot1) == -2*(ED(idxnot1 - 1) + DV(idxnot1 - 1) - EV(idxnot1 - 1) + M(idxnot1 - 1));

prob = optimproblem;
prob.Constraints.cons1 = cons1;
prob.Constraints.consd0 = consd0;
prob.Constraints.consED0 = consED0;
prob.Constraints.M0 = consM0;
prob.Constraints.S0 = consS0;
prob.Constraints.consED = consED;
prob.Constraints.consMnot1 = consMnot1;
prob.Constraints.consSnot1 = consSnot1;
prob.Constraints.consSend = S(end-1:end) == [0;0];
prob.Objective = 2*sum((DV(idxnot1) - DV(1:(N-1))).^2) + 1/2*sum((EV - DV).^2);

[sol,fval,eflag,output] = solve(prob)

Solving problem using lsqlin.

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.

sol = struct with fields:
    DV: [42×1 double]
    ED: [42×1 double]
     M: [42×1 double]
     S: [42×1 double]
     d: [42×1 double]

fval = 408.7581

eflag = 
    OptimalSolution

output = struct with fields:
            message: '↵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.↵↵<stopping criteria details>↵↵Optimization completed: The relative dual feasibility, 6.443259e-16,↵is less than options.OptimalityTolerance = 1.000000e-08, the complementarity measure,↵0.000000e+00, is less than options.OptimalityTolerance, and the relative maximum constraint↵violation, 1.268826e-16, is less than options.ConstraintTolerance = 1.000000e-08.↵↵'
          algorithm: 'interior-point'
      firstorderopt: 1.8041e-14
    constrviolation: 3.5527e-15
         iterations: 1
       linearsolver: 'sparse'
       cgiterations: []
             solver: 'lsqlin'

disp(sol.DV)

         0
    1.0536
    1.7606
    3.3011
    4.0630
    5.2401
    5.6296
    6.3319
    6.5257
    6.7626
    7.6048
    8.2659
    8.4143
    8.0901
    8.2154
    8.3246
    7.9480
    8.4945
    6.6038
    6.3063
    6.0307
    6.2112
    5.8960
    4.7592
    3.7700
    4.6838
    5.2324
    3.0558
    0.6131
   -1.2034
   -2.8448
   -4.2183
   -5.1642
   -3.9159
   -3.9084
   -5.3865
   -7.7170
   -8.9792
  -10.4856
  -12.1099
  -13.2549
  -14.2039

plot(1:N,EV','b-',1:N,sol.DV','r-')

You can tweak the objective function to give different results. I rather arbitrarily picked the multipliers 2 and 1/2 that you see in the objective definition.

Alan Weiss

MATLAB mathematical toolbox documentation

[Kevin Syc]

Hi Alan,

this solution looks very good! The only thing is, that I would need the DV to be composed of integers only. When I try to define the DV type to be integer, I get the following error:

Error using optim.problemdef.OptimizationProblem/solve
Quadratic problem with integer variables not supported.

Error in untitled (line 28)
[sol,fval,eflag,output] = solve(prob);

Is there a workaround for this or would this have to be done with a different optimizer?

[Alan Weiss]

Maybe you can call round on the answer?

Alan Weiss

MATLAB mathematical toolbox documentation

[Kevin Syc]

I could but then the array would no longer be true to satisfy the last two digits of the S to be 0 because the array that satisfies this condition is the one from the solver with decimal places.

[Kevin Syc]

Hi Alan, could this be solved with intlinprog to only give me integers for DV?

[Alan Weiss]

Sorry, at the moment I do not know of a satisfactory nonlinear integer solver that handles equality constraints. However, using some experimental code, I was able to obtain a reasonable-looking solution to your problem. This is sort of a magic trick, in that I am not free to show you the experimental code, but in case you are interested here are the results.

DV = [0     3     3     3     4     5     5     5     4     7     7     8     9    10     9     9     9     8   6     5     5     6     5     5     4     5     7     5     1    -3    -4    -7    -5    -3    -2    -5  -7   -10   -13   -13   -10   -10];
ED = [0     0     1     6     3     5     2     3     0    -2     1     0    -2    -3     1     2     1     6  -4     2     1     2     0    -4    -3     5     4    -5    -4    -1    -2    -2    -3     6     2    -2   -5     0    -2    -5    -4     0];
M = [0     0     1     7    10    15    17    20    20    18    19    19    17    14    15    17    18    24  20    22    23    25    25    21    18    23    27    22    18    17    15    13    10    16    18    16  11    11     9     4     0     0];
S = [0     0    -2   -14   -20   -30   -34   -40   -40   -36   -38   -38   -34   -28   -30   -34   -36   -48  -40   -44   -46   -50   -50   -42   -36   -46   -54   -44   -36   -34   -30   -26   -20   -32   -36   -32   -22   -22   -18    -8     0     0];
d = [0     1     5    -3     2    -3     1    -3    -2     3    -1    -2    -1     4     1    -1     5   -10  6    -1     1    -2    -4     1     8    -1    -9     1     3    -1     0    -1     9    -4    -4    -3  5    -2    -3     1     4     8];
EV = [0 2 -2 6 2 8 4 8 6 4 8 10 10 6 8 10 4 18 0 6 4 8 9 4 -4 6 16 4 -2 -2 -4 -6 -14 1 2 -2 -12 -8 -10 -14 -14 -18];
N = length(EV);
plot(1:N,EV,'b-',1:N,DV,'r-')

I hope to be able to show you the code in a few months. Sorry for the unsatisfactory response.

Alan Weiss

MATLAB mathematical toolbox documentation

[Kevin Syc]

Hi Alan thank you very much. I have one working code that goes through multiple permutations and it managed to find a good solution to this problem but doesn't find it for other problems. Problem is I can't do billions of permutations to see al of them and choose a good one. Thanks for your help though!

[Kevin Syc]

If I just might, add would you be able to give me a hand on setting up this problem in ga()? I think it might work for this problem but I can't make it work since I don't know how to input the constraints in.

[Alan Weiss]

Sorry, ga does not currently handle linear equality constraints. All of your equations are linear equality constraints.

Alan Weiss

MATLAB mathematical toolbox documentation

[Alan Weiss]

Now that R2021b has been released, I can show you how I was able to get the result I previously related as a "magic trick." The point is that in R2021b, ga accepts linear equality constraints.

EV = [0 2 -2 6 2 8 4 8 6 4 8 10 10 6 8 10 4 18 0 6 4 8 9 4 -4 6 16 4 -2 -2 -4 -6 -14 1 2 -2 -12 -8 -10 -14 -14 -18];
N = length(EV);
d = optimvar('d',1,N,"Type","integer","LowerBound",-20,"UpperBound",20);
ED = optimvar('ED',1,N,"Type","integer","LowerBound",-25,"UpperBound",25);
M = optimvar('M',1,N,"Type","integer","LowerBound",-40,"UpperBound",40);
S = optimvar('S',1,N,"Type","integer","LowerBound",-80,"UpperBound",20);
DV = optimvar('DV',1,N,"Type","integer","LowerBound",-50,"UpperBound",50);

prob = optimproblem;
idxnot1 = 2:N;
prob.Constraints.cons1 = d == DV - EV;
prob.Constraints.consd0 = d(1) == 0;
prob.Constraints.consED0 = ED(1) == 0;
prob.Constraints.M0 = M(1) == 0;
prob.Constraints.S0 = S(1) == 0;
prob.Constraints.consED = ED(idxnot1) == ED(idxnot1 - 1) + d(idxnot1 - 1);
prob.Constraints.consMnot1 = M(idxnot1) == ED(idxnot1) + M(idxnot1 - 1);
prob.Constraints.consSnot1 = S(idxnot1) == -2*(ED(idxnot1 - 1) + DV(idxnot1 - 1) - EV(idxnot1 - 1) + M(idxnot1 - 1));
prob.Constraints.consSend = S(end-1:end) == [0,0];
prob.Objective = 12*sum((DV(idxnot1) - DV(1:(N-1))).^2) + 1/2*sum((EV - DV).^2);

rng default
options = optimoptions("ga","PlotFcn","gaplotbestf","Display","none","MaxGenerations",100);
[sol,fval,eflag,output] = solve(prob,"Options",options);

disp(sol.DV)

Columns 1 through 11

0 3.0000 3.0000 3.0000 3.0000 5.0000 5.0000 5.0000 6.0000 7.0000 7.0000

Columns 12 through 22

8.0000 9.0000 9.0000 8.0000 9.0000 9.0000 8.0000 5.0000 6.0000 4.0000 6.0000

Columns 23 through 33

6.0000 6.0000 6.0000 6.0000 5.0000 5.0000 0 -2.0000 -5.0000 -6.0000 -4.0000

Columns 34 through 42

-3.0000 -3.0000 -5.0000 -8.0000 -11.0000 -14.0000 -12.0000 -9.0000 -8.0000

figure
plot(1:N,EV,'b-',1:N,sol.DV,'r-')

Alan Weiss

MATLAB mathematical toolbox documentation