Optimization with loops Error in optim.prob​lemdef.Opt​imizationP​roblem/sol​ve

15 views (last 30 days)
Q1. I want to perform certain operation on two large column vectors and calculate a cumulative sum in a certain order without using loops. For simplicity assume that I have two column vectors:
A=[1;3;5;7] and B = [10;12;14;19]
and I want to find a vector C = [1;1*10+3;1*10+3; (1*10+3)*12+5; ((1*10+3)*12+5)*14+7 ] = [1;13;161;2261]
What will be the MATLAB code to calculate vector C in MATLAB without using loops?
Q2. I have provided my code below and I have used the helpful suggestion given in the answer to use a loop in an optimization problem. The use of a optimization expression C in the code below is intended to be similar in terms of relevant elements to add and multiply as in the example given above but it contains the optimization variables. I get the error messages shown below when I run the code. I have also tried to write operations on C manually as shown in the last line of the code but I get the same error messages at the last line of the code Solution = solve(ALM).
How can I correct the code?
ERROR MESSAGES
Index exceeds the number of array elements (0).
Error in optim.internal.problemdef.SubsasgnExpressionImpl/computeLinearCoefficients
Error in optim.internal.problemdef.ExpressionTree/extractLinearCoefficients
Error in optim.internal.problemdef.ExpressionForest/extractLinearCoefficients
Error in optim.problemdef.OptimizationExpression/extractLinearCoefficients
Error in optim.problemdef.OptimizationConstraint/extractLinearCoefficients
Error in optim.problemdef.OptimizationProblem/compileConstraints
Error in prob2structImpl
Error in optim.problemdef.OptimizationProblem/solve
prices = [99.74 91.22 98.71 103.75 97.15];
cashFlows = [4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 102.5 5 4;4 5 0 105 104;4 105 0 0 0; 104 0 0 0 0];
obligations = [5 7 7 6 8 7 20 0]'*1000;
nt=size(cashFlows,1)
nb=size(cashFlows,2)
Rates = [0.01; 0.015; 0.017;0.019;0.02;0.025;0.027;0.029];
EndTimes = (1:nt)';
Disc = rate2disc(1,Rates,EndTimes);
%Number of bonds available
nBonds = [10;100;20;30;5]
ALM = optimproblem;
bonds = optimvar('bonds',nb,'Type','integer','LowerBound',0,'UpperBound',nBonds);
ALM.ObjectiveSense = 'minimize';
ALM.Objective = prices*bonds;
%Define the constraint
B=[1.01;1.020024752;1.02101183;1.02502363;1.024009823;1.050370059;1.039082218;1.043109481];
A=-(cashFlows*bonds-obligations);
C = optimexpr(nt,1);
C(1) = A(1);
for k = 2:numel(A)
C(k) = C(k-1) * B(k) + A(k);
end
ALM.Constraints.Const1 = (C.*Disc)/53.844 <=0.05;
%Solve the problem
Solution = solve(ALM)
%Here I show the intended operation on Optimization Expression C (in case
%for loop above does not operate as intended but I get the same error messages
C=[A(1);A(1)*B(1)+A(2);(A(1)*B(1)+A(2))*B(2)+A(3);(A(1)*B(1)+A(2))*B(2)+A(3))*B(3)+A(4);...
((A(1)*B(1)+A(2))*B(2)+A(3))*B(3)+A(4))*B(4)+A(5);
(((A(1)*B(1)+A(2))*B(2)+A(3))*B(3)+A(4))*B(4)+A(5))*B(5)+A(6);...
((((A(1)*B(1)+A(2))*B(2)+A(3))*B(3)+A(4))*B(4)+A(5))*B(5)+A(6))*B(6)+A(7);...
(((((A(1)*B(1)+A(2))*B(2)+A(3))*B(3)+A(4))*B(4)+A(5))*B(5)+A(6))*B(6)+A(7))*B(7)+A(8)];
  14 Comments
L Smith
L Smith on 10 Oct 2019
Edited: L Smith on 10 Oct 2019
B(i) in most cases (but not always) be greater than 1. 2^100 seems quite a big value
Matt J
Matt J on 11 Oct 2019
Edited: Matt J on 11 Oct 2019
That's if the B(i) are all less than 1.1, asin you rexample, it should not be a big issue.

Sign in to comment.

Accepted Answer

Matt J
Matt J on 10 Oct 2019
Edited: Matt J on 10 Oct 2019
I believe it will probably involve multiplication by some upper triangle matrices of 1 with some other operations but I cannot work it out.
Here it is,
d=cumprod(B(1:end-1));
d=[1;d(:)];
M=tril(d./d.');
ALM.Constraints.Const1 = ((M*A).*Disc)/53.844 <=0.05;
Solution = solve(ALM)
  2 Comments
Matt J
Matt J on 10 Oct 2019
Edited: Matt J on 10 Oct 2019
I don't know why the solver has trouble when a for-loop is used, but it isn't strictly related to the for-loop. For example, this version runs without error messages,
B=[1.01;1.020024752;1.02101183;1.02502363;1.024009823;1.050370059;1.039082218;1.043109481];
A=-(cashFlows*bonds-obligations);
C = optimexpr(nt,1);
C(1) = A(1);
for k = 2:numel(A)
C(k) = C(k-1) * B(k) + A(k);
end
ALM.Constraints.Const1 = C <= 5e6;
Solution = solve(ALM)
L Smith
L Smith on 10 Oct 2019
Thank you Matt.
I think this should be posted as a separate question as it appears to me an issue with the problem-based approach. Another reason why solver-based approach is better.

Sign in to comment.

More Answers (1)

Stephen23
Stephen23 on 8 Oct 2019
Edited: Stephen23 on 9 Oct 2019
I don't see a simple vectorized solution, but one for loop will be quite efficient:
A = [1;3;5;7]
B = [10;12;14;19]
C = nan(size(A)); % preallocate
C(1) = A(1);
for k = 2:numel(A)
C(k) = C(k-1) * B(k-1) + A(k);
end
Giving:
C =
1
13
161
2261
  5 Comments
Daniel M
Daniel M on 9 Oct 2019
I have a hunch there is a very clever solution using something like conv2 or filter, but it would be much more complicated to understand than the loop.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!