Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

State Estimation modifications for WLS #2414

Draft
wants to merge 6 commits into
base: develop
Choose a base branch
from
Draft

State Estimation modifications for WLS #2414

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
8b19466
Fix pwr injection results for merged buses in SE
marcopau Oct 7, 2024
0bbd43d
Added shunt results in SE
marcopau Oct 7, 2024
26ac3c7
Fix initialization bugs in SE
marcopau Oct 7, 2024
abc4619
Enhancements and logs to WLS estimator
marcopau Oct 7, 2024
e4b6fbf
FIx bugs with current mag meas in SE
marcopau Oct 7, 2024
d4342f0
Merge pull request #7 from marcopau/50Hz_developments
marcopau Oct 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion pandapower/auxiliary.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1746,7 +1746,7 @@ def _init_runse_options(net, v_start, delta_start, calculate_voltage_angles,
net._options = {}
_add_ppc_options(net, calculate_voltage_angles=calculate_voltage_angles,
trafo_model=trafo_model, check_connectivity=check_connectivity,
mode="pf", switch_rx_ratio=switch_rx_ratio, init_vm_pu=v_start,
mode="se", switch_rx_ratio=switch_rx_ratio, init_vm_pu=v_start,
init_va_degree=delta_start, enforce_q_lims=False, recycle=None,
voltage_depend_loads=False, trafo3w_losses=trafo3w_losses)
_add_pf_options(net, tolerance_mva="1e-8", trafo_loading="power",
Expand Down
111 changes: 83 additions & 28 deletions pandapower/build_bus.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,27 +42,29 @@ def ds_find(ar, bus): # pragma: no cover


@jit(nopython=True, cache=False)
def ds_union(ar, bus1, bus2, bus_is_pv): # pragma: no cover
def ds_union(ar, bus1, bus2, bus_is_pv, bus_is_active, merged_bus): # pragma: no cover
root1 = ds_find(ar, bus1)
root2 = ds_find(ar, bus2)
if root1 == root2:
return
if bus_is_pv[root2]:
if bus_is_active[root2] and ~bus_is_pv[root1]: # if root2 is an active bus (load or gen) and root1 is not a PV bus, it will merge root1 into root2
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here we should actually also make sure that root2 is not external grid or reference bus as well

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If root2 is an external grid it is still fine to merge root1 into root2. You want to keep the PV bus, and in this case root2 (which is the PV bus) would be kept.

ar[root1] = root2
merged_bus[root1] = True
else:
ar[root2] = root1
merged_bus[root2] = True


@jit(nopython=True, cache=False)
def ds_create(ar, switch_bus, switch_elm, switch_et_bus, switch_closed, switch_z_ohm,
bus_is_pv, bus_in_service): # pragma: no cover
bus_is_pv, bus_is_active, bus_in_service, merged_bus): # pragma: no cover
for i in range(len(switch_bus)):
if not switch_closed[i] or not switch_et_bus[i] or switch_z_ohm[i] > 0:
continue
bus1 = switch_bus[i]
bus2 = switch_elm[i]
if bus_in_service[bus1] and bus_in_service[bus2]:
ds_union(ar, bus1, bus2, bus_is_pv)
ds_union(ar, bus1, bus2, bus_is_pv, bus_is_active, merged_bus)


@jit(nopython=True, cache=False)
Expand All @@ -74,7 +76,7 @@ def fill_bus_lookup(ar, bus_lookup, bus_index):
bus_lookup[b] = bus_lookup[ar[ds]]


def create_bus_lookup_numba(net, bus_index, bus_is_idx, gen_is_mask, eg_is_mask):
def create_bus_lookup_numba(net, bus_index, bus_is_idx):
max_bus_idx = np.max(bus_index)
# extract numpy arrays of switch table data
switch = net["switch"]
Expand All @@ -86,18 +88,39 @@ def create_bus_lookup_numba(net, bus_index, bus_is_idx, gen_is_mask, eg_is_mask)
# create array for fast checking if a bus is in_service
bus_in_service = np.zeros(max_bus_idx + 1, dtype=bool)
bus_in_service[bus_is_idx] = True
# create mask arrays to distinguish different grid elements
_is_elements = net["_is_elements"]
eg_is_mask = _is_elements['ext_grid']
gen_is_mask = _is_elements['gen']
load_is_mask = _is_elements['load']
motor_is_mask = _is_elements['motor']
sgen_is_mask = _is_elements['sgen']
shunt_is_mask = _is_elements['shunt']
ward_is_mask = _is_elements['ward']
xward_is_mask = _is_elements['xward']
# create array for fast checking if a bus is pv bus
bus_is_pv = np.zeros(max_bus_idx + 1, dtype=bool)
bus_is_pv[net["ext_grid"]["bus"].values[eg_is_mask]] = True
bus_is_pv[net["gen"]["bus"].values[gen_is_mask]] = True
# create array for checking if a bus is active (i.e., it has some element connected to it)
bus_is_active = np.zeros(max_bus_idx + 1, dtype=bool)
bus_is_active[net["ext_grid"]["bus"].values[eg_is_mask]] = True
bus_is_active[net["gen"]["bus"].values[gen_is_mask]] = True
bus_is_active[net["load"]["bus"].values[load_is_mask]] = True
bus_is_active[net["motor"]["bus"].values[motor_is_mask]] = True
bus_is_active[net["sgen"]["bus"].values[sgen_is_mask]] = True
bus_is_active[net["shunt"]["bus"].values[shunt_is_mask]] = True
bus_is_active[net["ward"]["bus"].values[ward_is_mask]] = True
bus_is_active[net["xward"]["bus"].values[xward_is_mask]] = True
# create array that represents the disjoint set
ar = np.arange(max_bus_idx + 1)
merged_bus = np.zeros(len(ar), dtype=bool)
ds_create(ar, switch_bus, switch_elm, switch_et_bus, switch_closed, switch_z_ohm, bus_is_pv,
bus_in_service)
bus_is_active, bus_in_service, merged_bus)
# finally create and fill bus lookup
bus_lookup = -np.ones(max_bus_idx + 1, dtype=np.int64)
fill_bus_lookup(ar, bus_lookup, bus_index)
return bus_lookup
return bus_lookup, merged_bus


class DisjointSet(dict):
Expand Down Expand Up @@ -130,20 +153,35 @@ def create_consecutive_bus_lookup(net, bus_index):
return bus_lookup


def create_bus_lookup_numpy(net, bus_index, bus_is_idx, gen_is_mask, eg_is_mask,
closed_bb_switch_mask):
def create_bus_lookup_numpy(net, bus_index, closed_bb_switch_mask):
bus_lookup = create_consecutive_bus_lookup(net, bus_index)
net._fused_bb_switches = closed_bb_switch_mask & (net["switch"]["z_ohm"].values <= 0)
merged_bus = np.zeros(len(bus_lookup), dtype=bool)
if net._fused_bb_switches.any():
# Note: this might seem a little odd - first constructing a pp to ppc mapping without
# fused busses and then update the entries. The alternative (to construct the final
# mapping at once) would require to determine how to fuse busses and which busses
# are not part of any opened bus-bus switches first. It turns out, that the latter takes
# quite some time in the average usecase, where #busses >> #bus-bus switches.

# create mask arrays to distinguish different grid elements
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Redundant code

_is_elements = net["_is_elements"]
eg_is_mask = _is_elements['ext_grid']
gen_is_mask = _is_elements['gen']
load_is_mask = _is_elements['load']
motor_is_mask = _is_elements['motor']
sgen_is_mask = _is_elements['sgen']
shunt_is_mask = _is_elements['shunt']
ward_is_mask = _is_elements['ward']
xward_is_mask = _is_elements['xward']
# Find PV / Slack nodes -> their bus must be kept when fused with a PQ node
pv_list = [net["ext_grid"]["bus"].values[eg_is_mask], net["gen"]["bus"].values[gen_is_mask]]
pv_ref = np.unique(np.hstack(pv_list))
# Find active nodes -> their bus must be possibly kept when fused with a zero injection node
active_bus = [net["load"]["bus"].values[load_is_mask], net["motor"]["bus"].values[motor_is_mask], \
net["sgen"]["bus"].values[sgen_is_mask], net["shunt"]["bus"].values[shunt_is_mask], \
net["ward"]["bus"].values[ward_is_mask], net["xward"]["bus"].values[xward_is_mask]]
active_ref = np.unique(np.hstack(active_bus))
# get the pp-indices of the buses which are connected to a switch to be fused
fbus = net["switch"]["bus"].values[net._fused_bb_switches]
tbus = net["switch"]["element"].values[net._fused_bb_switches]
Expand All @@ -170,30 +208,48 @@ def create_bus_lookup_numpy(net, bus_index, bus_is_idx, gen_is_mask, eg_is_mask,
if any(i in fbus or i in tbus for i in pv_ref):
# for every disjoint set
for dj in disjoint_sets:
# check if pv buses are in the disjoint set dj
# check if pv and active buses are in the disjoint set dj
pv_buses_in_set = set(pv_ref) & dj
nr_pv_bus = len(pv_buses_in_set)
if nr_pv_bus == 0:
# no pv buses. Use any bus in dj
map_to = bus_lookup[dj.pop()]
else:
active_buses_in_set = set(active_ref) & dj
nr_active_bus = len(active_buses_in_set)
if nr_pv_bus > 0:
# one pv bus. Get bus from pv_buses_in_set
map_to = bus_lookup[pv_buses_in_set.pop()]
ref_bus = pv_buses_in_set.pop()
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not a very robust way since if pandapower decides to not have the first element as REF, this wont work and will be difficult to debug

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Explain what "if pandapower decides to not have the first element as REF" means. Feel free also to provide an example where this does not work, so I can look at it.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if the pv_buses_in_set first element which you are popping is not REF bus because of breaking change from pandapower

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In pv_buses_in_set you have the list of pv buses for the considered set. What I am popping is a PV bus. That bus is then taken as bus to be kept later when doing the merging among multiple buses.

else:
if nr_active_bus > 0:
# no pv bus but another active bus. Get bus from active_buses_in_set
ref_bus = active_buses_in_set.pop()
else:
# neither pv buses nor active buses. Use any bus in dj
ref_bus = dj.pop()
map_to = bus_lookup[ref_bus]
for bus in dj:
# update lookup
bus_lookup[bus] = map_to
if bus != ref_bus:
merged_bus[bus] = 1
else:
# no PV buses in set
for dj in disjoint_sets:
# use any bus in set
map_to = bus_lookup[dj.pop()]
active_buses_in_set = set(active_ref) & dj
nr_active_bus = len(active_buses_in_set)
if nr_active_bus > 0:
# no ov bus but another active bus. Get bus from active_buses_in_set
ref_bus = active_buses_in_set.pop()
else:
# neither pv buses nor active busese. Use any bus in dj
ref_bus = dj.pop()
map_to = bus_lookup[ref_bus]
for bus in dj:
# update bus lookup
bus_lookup[bus] = map_to
return bus_lookup
if bus != ref_bus:
merged_bus[bus] = 1
return bus_lookup, merged_bus


def create_bus_lookup(net, bus_index, bus_is_idx, gen_is_mask, eg_is_mask, numba):
def create_bus_lookup(net, bus_index, bus_is_idx, numba):
switches_with_pos_z_ohm = net["switch"]["z_ohm"].values > 0
if switches_with_pos_z_ohm.any() or not numba:
# if there are any closed bus-bus switches find them
Expand All @@ -208,12 +264,10 @@ def create_bus_lookup(net, bus_index, bus_is_idx, gen_is_mask, eg_is_mask, numba
net._impedance_bb_switches = np.zeros(switches_with_pos_z_ohm.shape)

if numba:
bus_lookup = create_bus_lookup_numba(net, bus_index, bus_is_idx,
gen_is_mask, eg_is_mask)
bus_lookup, merged_bus = create_bus_lookup_numba(net, bus_index, bus_is_idx)
else:
bus_lookup = create_bus_lookup_numpy(net, bus_index, bus_is_idx,
gen_is_mask, eg_is_mask, closed_bb_switch_mask)
return bus_lookup
bus_lookup, merged_bus = create_bus_lookup_numpy(net, bus_index, closed_bb_switch_mask)
return bus_lookup, merged_bus


def get_voltage_init_vector(net, init_v, mode, sequence=None):
Expand Down Expand Up @@ -301,11 +355,10 @@ def _build_bus_ppc(net, ppc, sequence=None):
bus_lookup = create_consecutive_bus_lookup(net, bus_index)
else:
_is_elements = net["_is_elements"]
eg_is_mask = _is_elements['ext_grid']
gen_is_mask = _is_elements['gen']
# eg_is_mask = _is_elements['ext_grid']
# gen_is_mask = _is_elements['gen']
bus_is_idx = _is_elements['bus_is_idx']
bus_lookup = create_bus_lookup(net, bus_index, bus_is_idx,
gen_is_mask, eg_is_mask, numba=numba)
bus_lookup, merged_bus = create_bus_lookup(net, bus_index, bus_is_idx, numba=numba)

n_bus_ppc = len(bus_index)
# init ppc with empty values
Expand Down Expand Up @@ -374,6 +427,8 @@ def _build_bus_ppc(net, ppc, sequence=None):

net["_pd2ppc_lookups"]["bus"] = bus_lookup
net["_pd2ppc_lookups"]["aux"] = aux
if mode != "nx":
net["_pd2ppc_lookups"]["merged_bus"] = merged_bus


def _build_bus_dc_ppc(net, ppc):
Expand Down
9 changes: 4 additions & 5 deletions pandapower/build_gen.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,9 +39,6 @@ def _build_gen_ppc(net, ppc):
mode = net["_options"]["mode"]
distributed_slack = net["_options"]["distributed_slack"]

if mode == "estimate":
return

_is_elements = net["_is_elements"]
gen_order = dict()
f = 0
Expand Down Expand Up @@ -211,6 +208,7 @@ def _build_pp_gen(net, ppc, f, t):
delta = net["_options"]["delta"]
gen_is = net._is_elements["gen"]
bus_lookup = net["_pd2ppc_lookups"]["bus"]
mode = net["_options"]["mode"]

gen_buses = bus_lookup[net["gen"]["bus"].values[gen_is]]
gen_is_vm = net["gen"]["vm_pu"].values[gen_is]
Expand All @@ -222,11 +220,12 @@ def _build_pp_gen(net, ppc, f, t):

# set bus values for generator buses
ppc["bus"][gen_buses[ppc["bus"][gen_buses, BUS_TYPE] != REF], BUS_TYPE] = PV
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why here all other buses are set to PV

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here there is no change. Gen buses are always modelled as PV nodes in pandapower, sgens as PQ.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes i am aware of that, the question was why here are the PQ buses exclusively excluded

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function only deals with gen_buses, and it assigns the label "PV" to them. Other buses are handled somewhere else.

ppc["bus"][gen_buses, VM] = gen_is_vm
if mode != "se":
ppc["bus"][gen_buses, VM] = gen_is_vm

add_q_constraints(net, "gen", gen_is, ppc, f, t, delta)
add_p_constraints(net, "gen", gen_is, ppc, f, t, delta)
if net._options["mode"] == "opf":
if mode == "opf":
# this considers the vm limits for gens
ppc = _check_gen_vm_limits(net, ppc, gen_buses, gen_is)
if "controllable" in net.gen.columns:
Expand Down
27 changes: 25 additions & 2 deletions pandapower/estimation/algorithm/base.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,9 @@ def __init__(self, tolerance, maximum_iterations, logger=std_logger):
self.r = None
self.H = None
self.hx = None
self.iterations = None
self.obj_func = None
logging.basicConfig(level=logging.DEBUG)

def estimate(self, eppci: ExtendedPPCI, **kwargs):
self.initialize(eppci)
Expand All @@ -83,6 +86,7 @@ def estimate(self, eppci: ExtendedPPCI, **kwargs):

current_error, cur_it = 100., 0
# invert covariance matrix
eppci.r_cov[eppci.r_cov<(10**(-6))] = 10**(-6)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This value of 10-6 should ideally be parameterizable

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is an option I am considering, but I am not yet sure if to put it as an optional parameter for the user or not (I see some pros and cons). In any case, this will be definitely improved in the next round of updates, for the moment it was urgent to avoid having 0 or close to 0 values in the covariance matrix.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is central to the feasibility of the solution and to evaluate how good the estimate is and where are problems in the network, so one could change it and then simulate multiple usecases

r_inv = csr_matrix(np.diagflat(1 / eppci.r_cov ** 2))
E = eppci.E
while current_error > self.tolerance and cur_it < self.max_iterations:
Expand All @@ -94,6 +98,14 @@ def estimate(self, eppci: ExtendedPPCI, **kwargs):
# jacobian matrix H
H = csr_matrix(sem.create_hx_jacobian(E))

# remove current magnitude measurements at the first iteration
# because with flat start they have null derivative
if cur_it == 0 and eppci.any_i_meas:
idx = eppci.idx_non_imeas
r_inv = r_inv[idx,:][:,idx]
r = r[idx,:]
H = H[idx,:]

# gain matrix G_m
# G_m = H^t * R^-1 * H
G_m = H.T * (r_inv * H)
Expand All @@ -106,17 +118,28 @@ def estimate(self, eppci: ExtendedPPCI, **kwargs):
E += d_E.ravel()
eppci.update_E(E)

# log data
current_error = np.max(np.abs(d_E))
obj_func = (r.T*r_inv*r)[0,0]
self.logger.debug("Current delta_x: {:.7f}".format(current_error))
self.logger.debug("Current objective function value: {:.1f}".format(obj_func))

# Restore full weighting matrix with current measurements
if cur_it == 0 and eppci.any_i_meas:
r_inv = csr_matrix(np.diagflat(1 / eppci.r_cov ** 2))

# prepare next iteration
cur_it += 1
current_error = np.max(np.abs(d_E))
self.logger.debug("Current error: {:.7f}".format(current_error))

except np.linalg.linalg.LinAlgError:
self.logger.error("A problem appeared while using the linear algebra methods."
"Check and change the measurement set.")
return False

# check if the estimation is successfull
self.check_result(current_error, cur_it)
self.iterations = cur_it
self.obj_func = obj_func
if self.successful:
# store variables required for chi^2 and r_N_max test:
self.R_inv = r_inv.toarray()
Expand Down
25 changes: 14 additions & 11 deletions pandapower/estimation/algorithm/matrix_base.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
from pandapower.pypower.idx_brch import F_BUS, T_BUS
from pandapower.pypower.dSbus_dV import dSbus_dV
from pandapower.pypower.dSbr_dV import dSbr_dV
from pandapower.pypower.dIbr_dV import dIbr_dV
from pandapower.pypower.dIbr_dV import dIbr_dV_new

from pandapower.estimation.ppc_conversion import ExtendedPPCI

Expand Down Expand Up @@ -150,16 +150,19 @@ def _dvabus_dV(V):

def _dimiabr_dV(self, V):
# for current we only interest in the magnitude at the moment
dif_dth, dif_dv, dit_dth, dit_dv, If, It = dIbr_dV(self.eppci['branch'], self.Yf, self.Yt, V)
dif_dth, dif_dv, dit_dth, dit_dv = map(lambda m: m.toarray(), (dif_dth, dif_dv, dit_dth, dit_dv))
difm_dth = (np.abs(1e-5 * dif_dth + If.reshape((-1, 1))) - np.abs(If.reshape((-1, 1))))/1e-5
difm_dv = (np.abs(1e-5 * dif_dv + If.reshape((-1, 1))) - np.abs(If.reshape((-1, 1))))/1e-5
ditm_dth = (np.abs(1e-5 * dit_dth + It.reshape((-1, 1))) - np.abs(It.reshape((-1, 1))))/1e-5
ditm_dv = (np.abs(1e-5 * dit_dv + It.reshape((-1, 1))) - np.abs(It.reshape((-1, 1))))/1e-5
difa_dth = (np.angle(1e-5 * dif_dth + If.reshape((-1, 1))) - np.angle(If.reshape((-1, 1))))/1e-5
difa_dv = (np.angle(1e-5 * dif_dv + If.reshape((-1, 1))) - np.angle(If.reshape((-1, 1))))/1e-5
dita_dth = (np.angle(1e-5 * dit_dth + It.reshape((-1, 1))) - np.angle(It.reshape((-1, 1))))/1e-5
dita_dv = (np.angle(1e-5 * dit_dv + It.reshape((-1, 1))) - np.angle(It.reshape((-1, 1))))/1e-5
difm_dth, difm_dv, ditm_dth, ditm_dv = dIbr_dV_new(self.eppci['branch'], self.Yf, self.Yt, V)
difm_dth, difm_dv, ditm_dth, ditm_dv = map(lambda m: m.toarray(), (difm_dth, difm_dv, ditm_dth, ditm_dv))
difa_dth, difa_dv, dita_dth, dita_dv = 0*difm_dth, 0*difm_dv, 0*ditm_dth, 0*ditm_dv
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why are all these values zero?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At the moment, we are not yet supporting current phase angles for state estimation. I will add a line of comment there.

# dif_dth, dif_dv, dit_dth, dit_dv, If, It = dIbr_dV(self.eppci['branch'], self.Yf, self.Yt, V)
# dif_dth, dif_dv, dit_dth, dit_dv = map(lambda m: m.toarray(), (dif_dth, dif_dv, dit_dth, dit_dv))
# difm_dth = (np.abs(1e-5 * dif_dth + If.reshape((-1, 1))) - np.abs(If.reshape((-1, 1))))/1e-5
# difm_dv = (np.abs(1e-5 * dif_dv + If.reshape((-1, 1))) - np.abs(If.reshape((-1, 1))))/1e-5
# ditm_dth = (np.abs(1e-5 * dit_dth + It.reshape((-1, 1))) - np.abs(It.reshape((-1, 1))))/1e-5
# ditm_dv = (np.abs(1e-5 * dit_dv + It.reshape((-1, 1))) - np.abs(It.reshape((-1, 1))))/1e-5
# difa_dth = (np.angle(1e-5 * dif_dth + If.reshape((-1, 1))) - np.angle(If.reshape((-1, 1))))/1e-5
# difa_dv = (np.angle(1e-5 * dif_dv + If.reshape((-1, 1))) - np.angle(If.reshape((-1, 1))))/1e-5
# dita_dth = (np.angle(1e-5 * dit_dth + It.reshape((-1, 1))) - np.angle(It.reshape((-1, 1))))/1e-5
# dita_dv = (np.angle(1e-5 * dit_dv + It.reshape((-1, 1))) - np.angle(It.reshape((-1, 1))))/1e-5
return difm_dth, difm_dv, ditm_dth, ditm_dv, difa_dth, difa_dv, dita_dth, dita_dv


Expand Down
Loading
Loading