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

8 views (last 30 days)
Kevin Syc on 29 Apr 2021
Edited: Alan Weiss on 23 Sep 2021 at 17:24
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.

Alan Weiss on 29 Apr 2021
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
Alan Weiss on 23 Sep 2021 at 15:15
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