3
$\begingroup$

The problem

I have a quadratically constrained quadratic optimisation problem (QCQP) which takes about 13 seconds to find the optimal solution. I want to iteratively optimise a large number of times (50-100 times) and clearly it will take over 10 minutes to run the same problem 50-100 times (with different initial conditions) but I feel that it can be made faster.

I have narrowed down the reason to FindMinimize being slow at solving QCQPs. Any suggestions on how can I speed things up? I have tried playing with AccuracyGoal and PrecisionGoal with values like 3 and 3 respectively in particular and they offer significant improvement (more than halved runtime) but since the optimal values are small numbers which later get multiplied by 100 when calculating other quantities, precision and accuracy both seem important and I would like to not compromise there if possible.

Miniature code with important part intact

I am presenting a very simple version of my code here with small optimisation horizon and other numbers since giving the entire code seems overkill (it is quite long) and one with larger optimisation horizons would require me to supply large ready-made expressions (which would otherwise be calculated by my actual code). The extremely long expressions are just data I have given so that you can run my code and so that you have an idea of how the constraints and objective function look like.

inputs1 = {u51, u52};
(*the objective function of optimisation problem*)
obj = 1.2326*10^-32 + (-0.6 + u51)^2 + (-0.6 + u52)^2 + (-0.6 + 
      100. u51) (3231.8 (-0.6 + 100. u51) + 
      4.28555 (0. - 753.982 u51 + 100 u52)) + (0. - 753.982 u51 + 
      100 u52) (4.28555 (-0.6 + 100. u51) + 
      3231.8 (0. - 753.982 u51 + 100 u52));

(*some constants to make sure you have all required values to run the
code*)
M = 2;
alpha = 0.2;
hor = 2;
P = {{6.78614, -0.962986, 4.91036, -0.839262, 
    1.84927, -0.218886}, {-0.962986, 1.8242, -1.07876, 
    2.54396, -0.536404, 1.83934}, {4.91036, -1.07876, 
    3.70497, -1.16611, 1.4347, -0.571865}, {-0.839262, 
    2.54396, -1.16611, 3.64699, -0.633915, 
    2.69184}, {1.84927, -0.536404, 1.4347, -0.633915, 
    0.56924, -0.362581}, {-0.218886, 1.83934, -0.571865, 
    2.69184, -0.362581, 2.02307}};

windowfunction = Array[HammingWindow, 2*M + 1, {-1/2, 1/2}];
ga = {{0.0071, -0.0564, 0.00580, 0.0104, 0.00240, -0.0321}}; 
de = {{0.0026}} ;


zetas = {{0., 0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0., 0.}, {0., 0., 
    0., 0., 0., 0.}, {0., 0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0., 
    0.}, {0., 0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0., 0.}, {0., 0., 
    0., 0., 0., 0.}, {0., 0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0., 
    0.}, {0., 0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0., 0.}, {0., 0., 
    0., 0., 0., 0.}, {0., 0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0., 
    0.}, {0., 0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0., 0.}, {0., 0., 
    0., 0., 0., 0.}, {0., 0., 0., 0., 0., 
    0.}, {0. + 0.0679348 (0. + 1. (0. + 100. u51)), 0., 0., 0., 0., 
    0.}, {0., 0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0., 0.}, {0., 0., 
    0., 0., 0., 0.}, {0. + 0.125 (0. + 1. (0. + 100. u51)), 0., 0., 
    0., 0., 0.}, {0. + 0.042875 (0. + 1. (0. + 100. u51)) + 
     0.0679348 (0. + 1. (0. - 753.982 u51 + 100. u52)), 
    0. + 0.25 (0. + 1. (0. + 100. u51)), 0., 0., 0., 0.}, {0., 0., 0.,
     0., 0., 0.}, {0., 0., 0., 0., 0., 
    0.}, {0. + 0.0679348 (0. + 1. (0. + 100. u51)), 0., 0., 0., 0., 
    0.}, {0. + 0.0233016 (0. + 1. (0. + 100. u51)) + 
     0.125 (0. + 1. (0. - 753.982 u51 + 100. u52)), 
    0. + 0.13587 (0. + 1. (0. + 100. u51)), 0., 0., 0., 
    0.}, {0. - 0.157768 (0. + 1. (0. + 100. u51)) + 
     0.0679348 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
     0.042875 (0. + 1. (0. - 753.982 u51 + 100. u52)), 
    0. + 0.0466033 (0. + 1. (0. + 100. u51)) + 
     0.25 (0. + 1. (0. - 753.982 u51 + 100. u52)), 
    0. + 0.13587 (0. + 1. (0. + 100. u51)), 0., 0., 0.}, {0, 0, 0, 0, 
    0, 0}, {0.0108696 (0. + 1. (0. + 100. u51)), 0., 0., 0., 0., 
    0.}, {0.00372826 (0. + 1. (0. + 100. u51)) + 
     0.0679348 (0. + 1. (0. - 753.982 u51 + 100. u52)), 
    0. + 0.0217391 (0. + 1. (0. + 100. u51)), 0., 0., 0., 
    0.}, {-0.0252429 (0. + 1. (0. + 100. u51)) + 
     0.125 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
     0.0233016 (0. + 1. (0. - 753.982 u51 + 100. u52)), 
    0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
     0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52)), 
    0. + 0.0217391 (0. + 1. (0. + 100. u51)), 0., 0., 
    0.}, {-0.0117118 (0. + 1. (0. + 100. u51)) + 
     0.042875 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
     0.0679348 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
     0.157768 (0. + 1. (0. - 753.982 u51 + 100. u52)), 
    0. - 0.0504859 (0. + 1. (0. + 100. u51)) + 
     0.25 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
     0.0466033 (0. + 1. (0. - 753.982 u51 + 100. u52)), 
    0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
     0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52)), 
    0. + 0.0217391 (0. + 1. (0. + 100. u51)), 0., 0.}};
outputs = {{0}, {0}, {0}, {0}, {0.}, {0.}, {0. + 
     1. (0. + 100. u51)}, {0. + 
     1. (0. - 753.982 u51 + 100. u52)}, {0. + 
     1. (0. + 0.98241 u51 - 0.132629 u52)}, {0. + 
     1. (0. + 0.131326 u51 - 0.0174146 u52)}, {0. + 
     1. (0. - 0.000346987 u51 + 0.0000464269 u52)}};

(*writing and aggregating constraints of the problem*)
domains = u51 \[Element] Reals && u52 \[Element] Reals;
setcons = {-15, -100} \[VectorLessEqual] {0., 
     0. + 100 u51} \[VectorLessEqual] {15, 100} && -10 <= u51 <= 
    10 && -10 <= u52 <= 10;
terminalcons = (-0.6 + 100. u51) (-0.00308956 + 
       1.35067*10^6 (-0.6 + 100. u51) - 
       2.39084*10^6 (0. - 753.982 u51 + 100 u52)) + 
    1.03986 (0.0224865 - 0.00297113 (-0.6 + 100. u51) + 
       0.00227077 (0. - 753.982 u51 + 100 u52)) + (0. - 753.982 u51 + 
       100 u52) (0.00236128 - 2.39084*10^6 (-0.6 + 100. u51) + 
       5.59947*10^6 (0. - 753.982 u51 + 100 u52)) <= 1;
(*next is the quadratic constraint calculation, which involves the
suspect Do loop*)
speccons = 


    0. + (0. + 1. (0. + 100. u51)) (0. + 
       5.11153*10^-8 (0. + 1. (0. + 100. u51))) <= 0.2 && 
  0. + (0. + 1. (0. + 100. u51)) (0. + 
       1.99669*10^-6 (0. + 1. (0. + 100. u51))) + (0. + 
       0.0679348 (0. + 1. (0. + 100. u51))) (0. + 
       6.78614 (0. + 0.0679348 (0. + 1. (0. + 100. u51)))) + (0. + 
       1. (0. - 753.982 u51 + 100. u52)) (0. + 
       1.60522*10^-6 (0. + 0.0679348 (0. + 1. (0. + 100. u51))) + 
       5.11153*10^-8 (0. + 1. (0. - 753.982 u51 + 100. u52))) + (0. + 
       0.0679348 (0. + 1. (0. + 100. u51))) (0. + 
       0.00005041 (0. + 0.0679348 (0. + 1. (0. + 100. u51))) + 
       1.60522*10^-6 (0. + 1. (0. - 753.982 u51 + 100. u52))) <= 0.2 &&
   0. + (0. + 1. (0. + 100. u51)) (0. + 
       6.76*10^-6 (0. + 1. (0. + 100. u51))) + (0. + 
       1. (0. - 753.982 u51 + 100. u52)) (0. + 
       0.0000100326 (0. + 0.125 (0. + 1. (0. + 100. u51))) + 
       1.99669*10^-6 (0. + 1. (0. - 753.982 u51 + 100. u52))) + (0. + 
       0.125 (0. + 1. (0. + 100. u51))) (0. + 
       0.00005041 (0. + 0.125 (0. + 1. (0. + 100. u51))) + 
       0.0000100326 (0. + 1. (0. - 753.982 u51 + 100. u52))) + (0. + 
       0.25 (0. + 1. (0. + 100. u51))) (0. + 
       1.8242 (0. + 0.25 (0. + 1. (0. + 100. u51))) - 
       0.962986 (0. + 0.042875 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.25 (0. + 1. (0. + 100. u51))) (0. + 
       0.00318096 (0. + 0.25 (0. + 1. (0. + 100. u51))) - 
       0.0000127513 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) - 
       0.00040044 (0. + 0.042875 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       1. (0. + 0.98241 u51 - 0.132629 u52)) (0. - 
       0.0000127513 (0. + 0.25 (0. + 1. (0. + 100. u51))) + 
       5.11153*10^-8 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       1.60522*10^-6 (0. + 0.042875 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.042875 (0. + 1. (0. + 100. u51)) + 
       0.0679348 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. - 
       0.00040044 (0. + 0.25 (0. + 1. (0. + 100. u51))) + 
       1.60522*10^-6 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       0.00005041 (0. + 0.042875 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.042875 (0. + 1. (0. + 100. u51)) + 
       0.0679348 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. - 
       0.962986 (0. + 0.25 (0. + 1. (0. + 100. u51))) + 
       6.78614 (0. + 0.042875 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. - 753.982 u51 + 100. u52)))) <= 0.2 &&
   0. + (0. + 1. (0. + 100. u51)) (0. + 
       1.99669*10^-6 (0. + 1. (0. + 100. u51))) + (0. + 
       1. (0. - 753.982 u51 + 100. u52)) (0. + 
       0.00001846 (0. + 0.0679348 (0. + 1. (0. + 100. u51))) + 
       6.76*10^-6 (0. + 1. (0. - 753.982 u51 + 100. u52))) + (0. + 
       0.0679348 (0. + 1. (0. + 100. u51))) (0. + 
       0.00005041 (0. + 0.0679348 (0. + 1. (0. + 100. u51))) + 
       0.00001846 (0. + 1. (0. - 753.982 u51 + 100. u52))) + (0. + 
       0.13587 (0. + 1. (0. + 100. u51))) (0. + 
       0.00318096 (0. + 0.13587 (0. + 1. (0. + 100. u51))) - 
       0.0000796957 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) - 
       0.00040044 (0. + 0.0233016 (0. + 1. (0. + 100. u51)) + 
          0.125 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       1. (0. + 0.98241 u51 - 0.132629 u52)) (0. - 
       0.0000796957 (0. + 0.13587 (0. + 1. (0. + 100. u51))) + 
       1.99669*10^-6 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       0.0000100326 (0. + 0.0233016 (0. + 1. (0. + 100. u51)) + 
          0.125 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.0233016 (0. + 1. (0. + 100. u51)) + 
       0.125 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. - 
       0.00040044 (0. + 0.13587 (0. + 1. (0. + 100. u51))) + 
       0.0000100326 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       0.00005041 (0. + 0.0233016 (0. + 1. (0. + 100. u51)) + 
          0.125 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.13587 (0. + 1. (0. + 100. u51))) (0. + 
       3.70497 (0. + 0.13587 (0. + 1. (0. + 100. u51))) + 
       4.91036 (0. - 0.157768 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.042875 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       1.07876 (0. + 0.0466033 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. - 
       0.157768 (0. + 1. (0. + 100. u51)) + 
       0.0679348 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       0.042875 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. + 
       4.91036 (0. + 0.13587 (0. + 1. (0. + 100. u51))) + 
       6.78614 (0. - 0.157768 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.042875 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.962986 (0. + 0.0466033 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. - 
       0.157768 (0. + 1. (0. + 100. u51)) + 
       0.0679348 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       0.042875 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. + 
       0.00004118 (0. + 0.13587 (0. + 1. (0. + 100. u51))) + 
       1.60522*10^-6 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) + 
       0.00005041 (0. - 0.157768 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.042875 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.00040044 (0. + 0.0466033 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.13587 (0. + 1. (0. + 100. u51))) (0. + 
       0.00003364 (0. + 0.13587 (0. + 1. (0. + 100. u51))) + 
       1.3113*10^-6 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) + 
       0.00004118 (0. - 0.157768 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.042875 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.00032712 (0. + 0.0466033 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       1. (0. + 0.131326 u51 - 0.0174146 u52)) (0. + 
       1.3113*10^-6 (0. + 0.13587 (0. + 1. (0. + 100. u51))) + 
       5.11153*10^-8 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) + 
       1.60522*10^-6 (0. - 0.157768 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.042875 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.0000127513 (0. + 0.0466033 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.0466033 (0. + 1. (0. + 100. u51)) + 
       0.25 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. - 
       0.00032712 (0. + 0.13587 (0. + 1. (0. + 100. u51))) - 
       0.0000127513 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
       0.00040044 (0. - 0.157768 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.042875 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
       0.00318096 (0. + 0.0466033 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.0466033 (0. + 1. (0. + 100. u51)) + 
       0.25 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. - 
       1.07876 (0. + 0.13587 (0. + 1. (0. + 100. u51))) - 
       0.962986 (0. - 0.157768 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.042875 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
       1.8242 (0. + 0.0466033 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. - 753.982 u51 + 100. u52)))) <= 0.2 && 
  0. + (0. + 1. (0. + 100. u51)) (0. + 
       5.11153*10^-8 (0. + 1. (0. + 100. u51))) + (0. + 
       1. (0. - 753.982 u51 + 100. u52)) (0. + 
       1.0905*10^-7 (0. + 1. (0. + 100. u51)) + 
       1.99669*10^-6 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
    0.0108696 (0. + 1. (0. + 100. u51)) (0. + 
       5.47935*10^-7 (0. + 1. (0. + 100. u51)) + 
       0.0000100326 (0. + 1. (0. - 753.982 u51 + 100. u52))) + (0. + 
       0.0217391 (0. + 1. (0. + 100. u51))) (0. + 
       0.00318096 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) - 
       0.00014664 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) - 
       0.00040044 (0.00372826 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       1. (0. + 0.98241 u51 - 0.132629 u52)) (0. - 
       0.00014664 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) + 
       6.76*10^-6 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       0.00001846 (0.00372826 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 
             1. (0. - 753.982 u51 + 100. u52)))) + (0.00372826 (0. + 
          1. (0. + 100. u51)) + 
       0.0679348 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. - 
       0.00040044 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) + 
       0.00001846 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       0.00005041 (0.00372826 (0. + 1. (0. + 100. u51)) + 
          0.0679348 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.0217391 (0. + 1. (0. + 100. u51))) (0. + 
       3.64699 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) - 
       0.839262 (-0.0117118 (0. + 1. (0. + 100. u51)) + 
          0.042875 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0679348 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
          0.157768 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
       2.54396 (0. - 0.0504859 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0466033 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       1.16611 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. - 
       0.0504859 (0. + 1. (0. + 100. u51)) + 
       0.25 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       0.0466033 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. + 
       2.54396 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) - 
       0.962986 (-0.0117118 (0. + 1. (0. + 100. u51)) + 
          0.042875 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0679348 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
          0.157768 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
       1.8242 (0. - 0.0504859 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0466033 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       1.07876 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 
             1. (0. - 753.982 u51 + 100. u52)))) + (-0.0252429 (0. + 
          1. (0. + 100. u51)) + 
       0.125 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       0.0233016 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. + 
       0.00004118 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) + 
       0.0000100326 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) + 
       0.00005041 (-0.0252429 (0. + 1. (0. + 100. u51)) + 
          0.125 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0233016 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.00040044 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.0217391 (0. + 1. (0. + 100. u51))) (0. + 
       0.00003364 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) + 
       8.19565*10^-6 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) + 
       0.00004118 (-0.0252429 (0. + 1. (0. + 100. u51)) + 
          0.125 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0233016 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.00032712 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. - 
       0.0504859 (0. + 1. (0. + 100. u51)) + 
       0.25 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       0.0466033 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. - 
       0.00058656 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) - 
       0.0000127513 (0. + 
          1. (0. - 0.000346987 u51 + 0.0000464269 u52)) - 
       0.00040044 (-0.0117118 (0. + 1. (0. + 100. u51)) + 
          0.042875 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0679348 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
          0.157768 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
       0.00318096 (0. - 0.0504859 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0466033 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.00032712 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       1. (0. + 0.131326 u51 - 0.0174146 u52)) (0. + 
       8.19565*10^-6 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) + 
       1.99669*10^-6 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) + 
       0.0000100326 (-0.0252429 (0. + 1. (0. + 100. u51)) + 
          0.125 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0233016 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.0000796957 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       1. (0. - 0.000346987 u51 + 0.0000464269 u52)) (0. + 
       2.3513*10^-6 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) + 
       5.11153*10^-8 (0. + 
          1. (0. - 0.000346987 u51 + 0.0000464269 u52)) + 
       1.60522*10^-6 (-0.0117118 (0. + 1. (0. + 100. u51)) + 
          0.042875 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0679348 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
          0.157768 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.0000127513 (0. - 0.0504859 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0466033 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
       1.3113*10^-6 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.00745652 (0. + 1. (0. + 100. u51)) + 
       0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. + 
       0.00006032 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) + 
       1.3113*10^-6 (0. + 
          1. (0. - 0.000346987 u51 + 0.0000464269 u52)) + 
       0.00004118 (-0.0117118 (0. + 1. (0. + 100. u51)) + 
          0.042875 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0679348 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
          0.157768 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.00032712 (0. - 0.0504859 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0466033 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
       0.00003364 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 
             1. (0. - 753.982 u51 + 100. u52)))) + (-0.0117118 (0. + 
          1. (0. + 100. u51)) + 
       0.042875 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       0.0679348 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
       0.157768 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. + 
       0.00007384 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) + 
       1.60522*10^-6 (0. + 
          1. (0. - 0.000346987 u51 + 0.0000464269 u52)) + 
       0.00005041 (-0.0117118 (0. + 1. (0. + 100. u51)) + 
          0.042875 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0679348 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
          0.157768 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.00040044 (0. - 0.0504859 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0466033 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
       0.00004118 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.0217391 (0. + 1. (0. + 100. u51))) (0. + 
       0.00010816 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) + 
       2.3513*10^-6 (0. + 
          1. (0. - 0.000346987 u51 + 0.0000464269 u52)) + 
       0.00007384 (-0.0117118 (0. + 1. (0. + 100. u51)) + 
          0.042875 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0679348 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
          0.157768 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.00058656 (0. - 0.0504859 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0466033 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
       0.00006032 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.00745652 (0. + 1. (0. + 100. u51)) + 
       0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. - 
       0.00032712 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) - 
       0.0000796957 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
       0.00040044 (-0.0252429 (0. + 1. (0. + 100. u51)) + 
          0.125 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0233016 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
       0.00318096 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52)))) + (0. + 
       0.00745652 (0. + 1. (0. + 100. u51)) + 
       0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. - 
       1.16611 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) + 
       4.91036 (-0.0117118 (0. + 1. (0. + 100. u51)) + 
          0.042875 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0679348 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
          0.157768 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       1.07876 (0. - 0.0504859 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0466033 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
       3.70497 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 
             1. (0. - 753.982 u51 + 100. u52)))) + (-0.0117118 (0. + 
          1. (0. + 100. u51)) + 
       0.042875 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
       0.0679348 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
       0.157768 (0. + 1. (0. - 753.982 u51 + 100. u52))) (0. - 
       0.839262 (0. + 0.0217391 (0. + 1. (0. + 100. u51))) + 
       6.78614 (-0.0117118 (0. + 1. (0. + 100. u51)) + 
          0.042875 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0679348 (0. + 1. (0. + 0.131326 u51 - 0.0174146 u52)) - 
          0.157768 (0. + 1. (0. - 753.982 u51 + 100. u52))) - 
       0.962986 (0. - 0.0504859 (0. + 1. (0. + 100. u51)) + 
          0.25 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
          0.0466033 (0. + 1. (0. - 753.982 u51 + 100. u52))) + 
       4.91036 (0. + 0.00745652 (0. + 1. (0. + 100. u51)) + 
          0.13587 (0. + 1. (0. - 753.982 u51 + 100. u52)))) <= 0.2;

(*finally finding the optimal solution*)
constraints = domains && speccons && terminalcons && setcons;
FindMinimum[{obj, constraints}, inputs1, Method -> "InteriorPoint", 
 MaxIterations -> 10000]

I am quite tired of this problem to be honest and would be extremely grateful for any suggestions. I have kept MaxIterations as 10000 because for larger optimisation horizon, this was needed to ensure I get an optimal solution satisfying tolerances. It seems though that it doesn't work for this smaller optimisation horizon.

Update to the question

Update 1: Based on Thorimur's comment and an experiment motivated by it, I have realised that my previous suspect, a Do loop, is not the culprit and hence I have modified this question to ask about how can FindMinimize be made faster (or what are the alternatives of course).

Update 2: Thank you very much for the answers, they all help. I am curious to know if using ConicOptimization function or QuadraticOptimization be faster than FindMinimum? Any comments would be really helpful. I thought along this line because all the answers here give things which speed up the code, but for horizon = 50, things still take quite long and it is practically impossible to post a miniature code here for horizon = 50.

Update 3: This is most likely too much to ask from you but I have uploaded my codes on GitHub here and with this, you can change the variable hor to whatever value you wish. The problem I wish to solve is hor=50 and M=20 and the miniature code I have given in the question (hor=2, M=2) is also derived from the uploaded codes. To run the notebook, all you will have to do is change the directory in SetDirectory to wherever you store the two packages (.wl files). In the off-chance you are not aware of this, the .nb file on GitHub looks ugly but upon pasting it in Mathematica, it will interpret the text and display it the way one would expect.

$\endgroup$
7
  • $\begingroup$ Removing Method -> "InteriorPoint" seems to do the trick! Do you have it there for a reason? (If this works, I'll edit my answer.) $\endgroup$ – thorimur Mar 30 at 20:15
  • $\begingroup$ It does do wonders for horizon = 2 but for some of the larger ones like horizon = 10, the difference is not more than a second :(. Since I was working with larger horizons and the method didn't make a difference, I just left it there but yeah I now see that it has implications. $\endgroup$ – ModCon Mar 30 at 20:27
  • $\begingroup$ Ah, I see. And I'm guessing combining it with Simplify (which seems safe, unlike FullSimplify) doesn't help? It seemed to halve the time in the given example, but I'm not sure how that generalizes. $\endgroup$ – thorimur Mar 30 at 20:36
  • $\begingroup$ Yea it isn't helping, doesn't scale well I think. It deteriorates the runtime from about 12 seconds to over a minute for horizon = 10. $\endgroup$ – ModCon Mar 30 at 20:45
  • $\begingroup$ Hmm. Unfortunately I don't know how to help further with FindMinimum; all I can suggest is that since you said you had to do it multiple different times, it might be worth trying to do those in parallel with e.g. Parallelize or other parallel functionality (though mathematica's parallelization can be finicky, and might require some extra manual work, e.g. distribution of definitions to the right kernels, ensuring things don't conflict, etc.). $\endgroup$ – thorimur Mar 30 at 22:27
3
$\begingroup$

I think the problem is not neither the Do loop or the fact that the code is not compiled. (I think FindMinimum will try to compile if it can.) The problem seems to be rather that FindRoot is not good at handling VectorLessEqual. This is a quite new feature and so I was surprised that it existed. Rephrasing the inequality constraints "in the good ol' way" seemed to help:

setcons = And @@ Thread[{-15, -100} <= {0., 0. + 100 u51} <= {15, 100}] && -10 <= u51 <= 10 && -10 <= u52 <= 10
simpconstraints = Simplify[speccons && terminalcons && setcons]; 
sol = FindMinimum[{obj, simpconstraints}, inputs1, MaxIterations -> 10000]; // AbsoluteTiming // First

simpconstraints /. sol[[2]]

-100 <= 0. + 100 u51 <= 100 && -10 <= u51 <= 10 && -10 <= u52 <= 10

0.034588

True So this is almost certainly a bug. Please file a bug report at Wolfram Support. By the way, I work with version 12.0.

PS.:

Also specifying the real numbers as domains seems to have complicated things more than it should. IIRC, FindMinimum assumes real numbers by default.

$\endgroup$
13
  • $\begingroup$ This is amazing! I was also searching for how to get rid of Do loops when defining constraints and And@@ with Table does the work and speeds things up even further. For horizon = 10, the runtime is now 1.4 seconds down from 12 seconds. However, my end problem has horizon = 50 and that takes an eternity even now. Could you please help me with that? My original code is in the question linked to this question so perhaps you would be able to find some culprits there. I use version 12.0 as well. Thank you! PS: I will file a bug report. $\endgroup$ – ModCon Mar 31 at 16:16
  • $\begingroup$ "However, my end problem has horizon = 50 and that takes an eternity even now. Could you please help me with that?" Maybe. But in the code posted above, hor does not do anything. I guess it is supposed to be used somehow to generate speccons. So it would be great if you updated your code. $\endgroup$ – Henrik Schumacher Mar 31 at 16:24
  • $\begingroup$ Umm, my code uses functions from a package I made and the variables which use those functions are literally over 300 rows long so I am not sure if I should copy paste such large matrices here. Instead, would the code given here ( mathematica.stackexchange.com/q/243814/73262) work? I am hoping improvements made there for hor = 2 would be scalable, just like yours in this answer are. Sorry for the inconvenience. $\endgroup$ – ModCon Mar 31 at 17:13
  • $\begingroup$ It might be too much to ask but in case you would like to, I have given a GitHub link to where my codes are stored in 'Update 3' of my question. You can then try hor = 50 as well. $\endgroup$ – ModCon Apr 1 at 12:44
  • 1
    $\begingroup$ Ah, a github link is a very good idea. I am a bit busy at the moment though. Maybe I'll find some time on the long weekend ahead. $\endgroup$ – Henrik Schumacher Apr 1 at 12:53
2
$\begingroup$

Note: I found that this answer doesn't actually work, but I'll leave it here as it provides some insight nonetheless.

Upon surrounding different sections with Echo[AbsoluteTiming[...], "section-label", First][[2]], your Do loop actually seems quite fast (0.002 seconds on my machine). FindMinimum, however, takes quite a while. (To test this more easily, simply remove FindMinimum and see how long it takes.)

The definition of speccons is evaluated once, and then the output is used in FindMinimum, so it is not as though we are running the Do loop in the definition each time we iterate in FindMinimum. speccons, as a fully evaluated expression, is nonetheless used in the process of FindMinimum, though, and it is quite complicated, so we might ask if there's a way to simplify it. And indeed, modifying constraints to be

constraints = FullSimplify[domains && speccons && terminalcons && setcons];

seems to work wonders!

EDIT: No it doesn't! FullSimplify does not actually preserve truth values when manipulating logical formulae, apparently. It seems FullSimplify may actually cause a different minimum to be found. Indeed,

constraints = 
  FullSimplify[domains && speccons && terminalcons && setcons];
Simplify[domains && speccons && terminalcons && setcons /. 
  Last@FindMinimum[{obj, constraints}, inputs1, 
    Method -> "InteriorPoint", MaxIterations -> 10000]]

shows us that the minimum found with FullSimplify actually violates the constraints.

We can check that the resulting constraints are actually equivalent. Simplify seems to return unquestionably equivalent constraints, but no speedup, whereas FullSimplify's constraints might be weaker; To see the logical relationships, consider the outputs of

Simplify@
 Equivalent[Simplify[domains && speccons && terminalcons && setcons], 
  domains && speccons && terminalcons && setcons]

Simplify@
 Equivalent[FullSimplify[domains && speccons && terminalcons && setcons], 
  domains && speccons && terminalcons && setcons]

Simplify@
 Implies[domains && speccons && terminalcons && setcons, 
  FullSimplify[domains && speccons && terminalcons && setcons]]
$\endgroup$
1
  • $\begingroup$ Thanks for this piece of work! I will keep this idea handy for future. $\endgroup$ – ModCon Mar 30 at 20:34
2
$\begingroup$

This works in 0.25 seconds.

zetas=FullSimplify[zetas];
terminalcons=FullSimplify[terminalcons];
speccons=List@@{FullSimplify[speccons]/.(1.0->1)};
setcons={-1.0<=u51<=1.0,-10.0<=u52<=10.0};
FindMinimum[Flatten@{obj, speccons, terminalcons, setcons}, {u51, u52}]
(* {0.660595,{u51->0.00600015,u52->0.04524}} *)
$\endgroup$
7
  • $\begingroup$ Have you missed out on -100 <= 0. + 100 u51 <= 100 in setcons? Without it, it does take 0.25 seconds but with it, it takes around 34 seconds for me. $\endgroup$ – ModCon Apr 1 at 8:59
  • $\begingroup$ But these do help even then! $\endgroup$ – ModCon Apr 1 at 9:07
  • $\begingroup$ -100 <= 0. + 100 u51 <= 100 is equivalent to -1.0<=u51<=1.0 which I did include. $\endgroup$ – Ted Ersek Apr 2 at 0:16
  • $\begingroup$ Ah sorry, should have noticed that. Also, what is the reason for replacing 1.0 with 1? I am curious because point #1 here (blog.wolfram.com/2011/12/07/…) recommends using floating point numbers. $\endgroup$ – ModCon Apr 2 at 7:54
  • $\begingroup$ Replacing 1.0 with 1 changes 1.0*u51 into u51 which is the same thing and a bit more efficient. If u51 isn't approximate to start it will be converted to one soon after because of all the other approximate numbers. I could have also changes 0.0 to 0 in terminalcons. As it is terminalcons has (0.0 - 753.982*u51 + 100*u52) when it could be ( -753.982*u51 + 100*u52). Also, since floating point numbers are faster we could change (100*u52) into (100.0*u52). $\endgroup$ – Ted Ersek Apr 3 at 20:46
1
$\begingroup$

Try

setcons = And @@ Thread[{-15, -100} <= {0., 0. + 100 u51} <={15,100}] && -10 <= u51 <= 10 && -10 <= u52 <= 10
constraints = speccons && terminalcons && setcons;
R = ImplicitRegion[constraints, {u51, u52}];
NArgMin[obj, {u51, u52} \[Element] R]

(* {0.00600015, 0.0452401} *)

RegionPlot[R]

enter image description here

$\endgroup$
1
  • $\begingroup$ It works! But whenever I use inputs1={u51,u52} and then use inputs1 instead of {u51,u52} in ImplicitRegion and NArgMin, it does not solve the problem and I need to do this for large horizons. Any idea how to circumvent this? $\endgroup$ – ModCon Apr 1 at 10:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.