search

Electa-Git / PowerModelsDistributionStateEstimation.jl

A Julia Package for Power System State Estimation.
BSD 3-Clause "New" or "Revised" License
32 stars 12 forks source link

Discuss best numerical approach for arc tangent constraint #7

Closed MartaVanin closed 4 years ago

MartaVanin commented 4 years ago

Among the measurement conversion possibilities, we have some that involve arctangents.

E.g.: current_angle_branch = tan^-1(q/p) The following won't work, because nothing prevents p[c] from becoming 0:

    JuMP.@NLconstraint(pm.model, [c in 1:nph],
        current_angle_branch_variable[component_id][c] == atan(q[c]/p[c])
        )

Current workaround (which works):

    JuMP.@NLconstraint(pm.model, [c in 1:nph],
        current_angle_branch_variable[component_id][c] == atan(q[c]/p[c]+0.00001)
        )

Can we do something more elegant? Also, I am a bit reluctant about the atan(x) function: it returns a value in the interval [-pi/2, pi/2], while we do have angles in the second/third quadrant, so this approach might compromise the estimation.

There exists atan(y,x) in julia, equivalent to atan(z) above, where z = y/x and the returned value is in [-pi, pi], which is ideal. But it is not available in JuMP. It is not available in JuMP because JuMP imports derivative from Calculus.jl, and nobody has written the derivative for atan(x,y), while it exists for atan(x) (see Calculus.jl/src/differentitate.jl line 141). I suppose it is not immediate to generalize that derivative to add a + or - \pi when x is negative, so I see only one work-around: we make a guess on where the angles will be located. For example, being distrib. networks quite resistive, we assume that current (phasor) angles will be in the same quadrant as the voltage's (where voltage angles are around [0, 120, -120] degrees). Then, the constraint becomes:

    JuMP.@NLconstraint(pm.model, 
        current_angle_branch_variable[component_id][1] == atan(q[1]/p[1]+0.00001)
        )
    JuMP.@NLconstraint(pm.model, 
        current_angle_branch_variable[component_id][2] == atan(q[2]/p[2]+0.00001)+pi
        )
    JuMP.@NLconstraint(pm.model, 
        current_angle_branch_variable[component_id][3] == atan(q[3]/p[3]+0.00001)-pi
        )

This is prone to poor accuracy and convergence issues if something goes wrong, but we can say that if the user wants better than this, he should use a variable space that does not rely on an arctan conversion given his measured values.

timmyfaraday commented 4 years ago

isn't the ca = va + theta = va + tan-1(q,p)

timmyfaraday commented 4 years ago

http://www-labs.iro.umontreal.ca/~mignotte/IFT2425/Documents/EfficientApproximationArctgFunction.pdf

MartaVanin commented 4 years ago

isn't the ca = va + theta = va + tan-1(q,p) yeah, it's ca = va - tan-1(q,p), actually. The problem remains, but I can try registering a function in JuMP based on the link you sent. I will try to wrap up all the other transformations and will let you know.

timmyfaraday commented 4 years ago

Doesn't the sign depend on whether the current is 'inductive' or 'capacitive'? They are approximations, some seem to be very accurate... I suggest you pick one of them, get it working, and move on... Later on, we'll look into improving the formulation!

MartaVanin commented 4 years ago

closed in the merge. solved by "preprocessing" measurement conversion. not sure why the merge didn't close this automatically.