Sedumi Numerical error!!!!

[Muhammad Usman]

Dear Johan,

My solver (Sedumi) keeps falling into numerical error but I am unable to understand why. I want to clarify some concepts regarding the ranges of constraints and value of primal-dual residuals. I have few questions as reported below

Q1. The coefficient ranges of equality constraints are shown in red colour. Is this range large enough for the solver that it cannot solve the problem?

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|    ID|                          Constraint|     Coefficient range|               Tag|

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|    #1|             Equality constraint 1x1|         1 to 198.7109|      Active PQ L5|

|    #2|             Equality constraint 1x1|   0.48432 to 198.7109|    Reactive PQ L5|

|    #3|             Equality constraint 1x1|    25.098 to 198.7109|      Active PQ L6|

|    #4|             Equality constraint 1x1|    25.098 to 198.7109|    Reactive PQ L6|

|    #5|             Equality constraint 1x1|    25.098 to 198.7109|      Active PQ L7|

|    #6|             Equality constraint 1x1|    25.098 to 198.7109|    Reactive PQ L7|

|    #7|             Equality constraint 1x1|         1 to 198.7109|     Active Ntrl 8|

|    #8|             Equality constraint 1x1|   0.48432 to 198.7109|   Reactive Ntrl 8|

|    #9|         Element-wise inequality 2x1|      0.9025 to 1.1025|       Voltage L 5|

|   #10|         Element-wise inequality 2x1|      0.9025 to 1.1025|       Voltage L 6|

|   #11|         Element-wise inequality 2x1|      0.9025 to 1.1025|       Voltage L 7|

|   #12|   Equality constraint (complex) 4x4|              0.5 to 1|     Voltage Slack|

|   #13|     Matrix inequality (complex) 8x8|                1 to 1|                  |

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Q2. The objective function is being generated using a for loop. 

for k=1:4

objective = objective+(trace(Pnode{k,1}*X)+trace(Qnode{k,1}*X));

end

The individual coefficient ranges of the term  (trace(Pnode{k,1}*X)+trace(Qnode{k,1}*X)) are 

[k=1; Coeffiecient range: 10.8023 to 314.0324

k=2;  Coeffiecient range: 10.8023 to 314.0324

k=3;  Coeffiecient range: 10.8023 to 314.0324

k=4;  Coeffiecient range: 10.8023 to 314.0324]

However, after the loop is run, I got the following coefficient range of my objective function. 

Objective

Linear scalar (real, 48 variables)

Coeffiecient range: 3.5527e-15 to 314.0324

Following is the output after each iteration

Linear scalar (real, 15 variables)

Coeffiecient range: 10.8023 to 314.0324

Linear scalar (real, 28 variables)

Coeffiecient range: 7.1054e-15 to 314.0324

Linear scalar (real, 39 variables)

Coeffiecient range: 3.5527e-15 to 314.0324

Linear scalar (real, 48 variables)

Coeffiecient range: 3.5527e-15 to 314.0324

It can be seen that in the second run of the loop when the current value of the objective is added to the previous value of objective, the coefficient ranges get weird. Can you tell me why is this happening? and how can I avoid it? Maybe, this is the reason why am I not getting the result?

Q3: For inequality constraint, I know that primal residual value should not be negative. For equality constraint, the residual should be zero. Am I right? The last constraint (#13) is X>=0 (i.e, X should be PSD). Is it correct that if a system is solved succesfully, primal residual for such constraint would be zero?

 

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|    ID|               Constraint|   Primal residual|   Dual residual|               Tag|

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|    #1|      Equality constraint|       -9.1129e-09|             NaN|      Active PQ L5|

|    #2|      Equality constraint|       -4.4135e-09|             NaN|    Reactive PQ L5|

|    #3|      Equality constraint|       -7.9851e-14|             NaN|      Active PQ L6|

|    #4|      Equality constraint|        -3.793e-14|             NaN|    Reactive PQ L6|

|    #5|      Equality constraint|       -3.2863e-14|             NaN|      Active PQ L7|

|    #6|      Equality constraint|       -9.3286e-14|             NaN|    Reactive PQ L7|

|    #7|      Equality constraint|       -9.1129e-09|             NaN|     Active Ntrl 8|

|    #8|      Equality constraint|       -4.4136e-09|             NaN|   Reactive Ntrl 8|

|    #9|   Elementwise inequality|          0.073952|             NaN|       Voltage L 5|

|   #10|   Elementwise inequality|          0.098821|             NaN|       Voltage L 6|

|   #11|   Elementwise inequality|          0.091293|             NaN|       Voltage L 7|

|   #12|      Equality constraint|       -2.4715e-08|             NaN|     Voltage Slack|

|   #13|        Matrix inequality|       -1.4586e-08|             NaN|                  |

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Regards

Muhammad

[Muhammad Usman]

Hi Johan, 

I have figured out why I have numerical garbage in my objective function. Since, it is the summation of Pnode(1,1)+Pnode(2,1)+Pnode(3,1)+....In this summation, few terms, when added, should give zero but MATLAB is not giving hard zero. Rather, it is producing 7.1054e-15 value. That's why I have this coefficient range. Can you help me how to solve this problem? 

[Johan Löfberg]

Both coefficient range and residuals look fairly normal. Solvers work with finite precision, and better than 1e-8 is typically not to be expected

SeDuMi is just being honest about its numerical progress. You should try Mosek, it is typically more numerically robust.

[Muhammad Usman]

Hi Johan, thank you for the reply. Any suggestions about numerical garbage value in the objective function?

[Johan Löfberg]

It depends on the context. For instance, if the model is 1e-15*x+(y+z)^2 s.t 0*x+y+z==1 (and that 1e-15 obviously was supposed to be 0), you should clean it (using clean if you cannot do it manually), but if it is 1e-15*x^2+(y+z)^2 s.t 0*x+y+z==1 it clearly makes no difference really , etc.