4

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.3506710^6 (-0.6 + 100. u51) - 
       2.3908410^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.3908410^6 (-0.6 + 100. u51) + 
       5.5994710^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))) &lt;= 0.2 &amp;&amp; 



    

    1. 

    • (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. +
        2. 

      
1.6052210^-6 (0. + 0.0679348 (0. + 1. (0. + 100. u51))) + 
5.1115310^-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 &&
      • 

    

    2. 

    2. 

    • (0. + 1. (0. + 100. u51)) (0. +
      • 

    
6.76*10^-6 (0. + 1. (0. + 100. u51))) + (0. +
      

    1. (0. - 753.982 u51 + 100. u52)) (0. +
      2. 

    
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. -
      2. 

    
0.0000127513 (0. + 0.25 (0. + 1. (0. + 100. u51))) + 
 5.1115310^-8 (0. + 1. (0. + 0.98241 u51 - 0.132629 u52)) + 
 1.6052210^-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 &&
    3. 

    3. 

    • (0. + 1. (0. + 100. u51)) (0. +
      • 

    
1.99669*10^-6 (0. + 1. (0. + 100. u51))) + (0. +
      

    1. (0. - 753.982 u51 + 100. u52)) (0. +
      2. 

    
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. -
      2. 

    
0.0000796957 (0. + 0.13587 (0. + 1. (0. + 100. u51))) + 
 1.9966910^-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.6052210^-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. +
      2. 

    
1.311310^-6 (0. + 0.13587 (0. + 1. (0. + 100. u51))) + 
 5.1115310^-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 &&
    4. 

    4. 

    • (0. + 1. (0. + 100. u51)) (0. + 
5.11153*10^-8 (0. + 1. (0. + 100. u51))) + (0. +
        

      1. (0. - 753.982 u51 + 100. u52)) (0. +
        2. 

      
1.090510^-7 (0. + 1. (0. + 100. u51)) + 
1.9966910^-6 (0. + 1. (0. - 753.982 u51 + 100. u52))) +
      • 

    

    5. 



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)))) &lt;= 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.

ModCon
  • 245
  • 1
  • 5
  • Removing Method -> "InteriorPoint" seems to do the trick! Do you have it there for a reason? (If this works, I'll edit my answer.) – thorimur Mar 30 '21 at 20:15
  • 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. – ModCon Mar 30 '21 at 20:27
  • 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. – thorimur Mar 30 '21 at 20:36
  • 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. – ModCon Mar 30 '21 at 20:45
  • 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.). – thorimur Mar 30 '21 at 22:27
  • Making use of NMinimize[{obj, constraints}, inputs1, Method -> "DifferentialEvolution", MaxIterations -> 10000], I obtain error communications. – user64494 Apr 01 '21 at 14:27
  • @user64494 It does give an error about not being able to find any feasible initial point for me as well, but the code doesn't stop running there and eventually gives an answer very close to one by FindMinimum. Any particular reason you chose this method? I avoided NMinimize because my problem is convex so I can use the faster FindMinimum instead and that would give me global minimizer as well. – ModCon Apr 01 '21 at 15:11

4 Answers4

4

I think the problem is 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.

Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309
  • 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.

     – ModCon Mar 31 '21 at 16:16
  • "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. – Henrik Schumacher Mar 31 '21 at 16:24
  • 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 ( https://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. – ModCon Mar 31 '21 at 17:13
  • 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. – ModCon Apr 01 '21 at 12:44
  • 1
    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. – Henrik Schumacher Apr 01 '21 at 12:53
  • Okay, I had a look into the git repo. Thank you for creating it. The code there is indeed far more complex than I can understand it in the bit of time I can afford to invest. The issue there is that you try to create lists of constraints that might be very redundant while the "InteriorPoint" method employs all of them. – Henrik Schumacher Apr 02 '21 at 08:41
  • Moreover and more importantly, you try to handle the constraints in a symbolic way. That is fine for a dozen constraints, but it seems to go awfully wrong here. For example, the expression speccons that you generate in the function optimalinputs1 has a LeafCount of 214598213. That's humongous. – Henrik Schumacher Apr 02 '21 at 08:41
  • If I had to solve these problems, then I would rephrase the problem into functional form: Instead of an expression for the objective, I would write an efficient function that eats a numerical vector and returns the function value. I would do the same with the derivative. And the same with the equality constraints and the inequality constraints and their derivatives separately. Then I would set up the penalty/barrier formulation myself (see https://reference.wolfram.com/language/tutorial/ConstrainedOptimizationLocalNumerical.html#442312913). – Henrik Schumacher Apr 02 '21 at 08:54
  • Once the penalty/barrier formulation is set up (with an efficient implemenation of the objective function and its gradient), an unconstrained optimization method of FindMinimum can be employed, an the gradient can be specified with the Gradient option of FindMinimum. Unfortunately, FindMinimum does not seem to support a corresponding way to communicate the derivatives of the constraints... – Henrik Schumacher Apr 02 '21 at 08:55
  • Sorry, I had some health issues hence the delay in responding. This is very interesting and I will try this out. I am also curious if using a dedicated function like ConicOptimisation would be faster than using the unconstrained FindMinimum as you have said. Could you please share your thoughts on this? – ModCon Apr 04 '21 at 23:03
  • Oh, I hope you feel well again! Hm. I do not have any experience with ˋConicOptimizationˋ or any of the other new dedicated optimization tools (ˋQuadraticOptimizationˋ, ˋSemidefiniteOptimizationˋ, etc.) If you. an reformulated your problem into one of these, then they will probably work much better. However, IIRC, they work only with linear objective functions... =/ – Henrik Schumacher Apr 05 '21 at 06:00
  • Thank you, I do feel well again. This sounds promising, I will give it a try by probably converting the quadratic problem to a linear one using slack variables. – ModCon Apr 05 '21 at 15:05
  • Ah, that sounds like a good idea! – Henrik Schumacher Apr 05 '21 at 15:06
3

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]]
thorimur
  • 9,010
  • 18
  • 32
3

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}} *)
Ted Ersek
  • 7,124
  • 19
  • 41
  • 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. – ModCon Apr 01 '21 at 08:59
  • But these do help even then! – ModCon Apr 01 '21 at 09:07
  • -100 <= 0. + 100 u51 <= 100 is equivalent to -1.0<=u51<=1.0 which I did include. – Ted Ersek Apr 02 '21 at 00:16
  • Ah sorry, should have noticed that. Also, what is the reason for replacing 1.0 with 1? I am curious because point #1 here (https://blog.wolfram.com/2011/12/07/10-tips-for-writing-fast-mathematica-code/) recommends using floating point numbers. – ModCon Apr 02 '21 at 07:54
  • Replacing 1.0 with 1 changes 1.0u51 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.982u51 + 100u52) when it could be ( -753.982u51 + 100u52). Also, since floating point numbers are faster we could change (100u52) into (100.0*u52). – Ted Ersek Apr 03 '21 at 20:46
  • FindMinimum uses algebra and stuff like 0.0+1.0*u51 only adds extra steps that don't change the outcome. Replace 0.0 with 0 and 1.0 with 1 and that is simply u51. – Ted Ersek Apr 04 '21 at 12:50
  • Ah, I get it now. Thanks! – ModCon Apr 04 '21 at 23:16
1

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

Cesareo
  • 3,963
  • 7
  • 11
  • 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? – ModCon Apr 01 '21 at 10:07