Evaluates Hessian of Lagrangian for AC OPF.

def _hessfcn(self, x, lmbda):
        """ Evaluates Hessian of Lagrangian for AC OPF.
        """
        Pgen = x[self._Pg.i1:self._Pg.iN + 1] # Active generation in p.u.
        Qgen = x[self._Qg.i1:self._Qg.iN + 1] # Reactive generation in p.u.

        for i, g in enumerate(self._gn):
            g.p = Pgen[i] * self._base_mva # active generation in MW
            g.q = Qgen[i] * self._base_mva # reactive generation in MVAr

        Vang = x[self._Va.i1:self._Va.iN + 1]
        Vmag = x[self._Vm.i1:self._Vm.iN + 1]
        V = Vmag * exp(1j * Vang)
        nxtra = self._nxyz - 2 * self._nb

        #------------------------------------------------------------------
        #  Evaluate d2f.
        #------------------------------------------------------------------

        d2f = self._d2f(x) * self.opt["cost_mult"]
        # TODO: Generalised cost model.

        #------------------------------------------------------------------
        #  Evaluate Hessian of power balance constraints.
        #------------------------------------------------------------------

        nlam = len(lmbda["eqnonlin"]) / 2
        lamP = lmbda["eqnonlin"][:nlam]
        lamQ = lmbda["eqnonlin"][nlam:nlam + nlam]
        Gpaa, Gpav, Gpva, Gpvv = self.om.case.d2Sbus_dV2(self._Ybus, V, lamP)
        Gqaa, Gqav, Gqva, Gqvv = self.om.case.d2Sbus_dV2(self._Ybus, V, lamQ)

        d2G = vstack([
            hstack([
                vstack([hstack([Gpaa, Gpav]),
                        hstack([Gpva, Gpvv])]).real +
                vstack([hstack([Gqaa, Gqav]),
                        hstack([Gqva, Gqvv])]).imag,
                csr_matrix((2 * self._nb, nxtra))]),
            hstack([
                csr_matrix((nxtra, 2 * self._nb)),
                csr_matrix((nxtra, nxtra))
            ])
        ], "csr")

        #------------------------------------------------------------------
        #  Evaluate Hessian of flow constraints.
        #------------------------------------------------------------------

        nmu = len(lmbda["ineqnonlin"]) / 2
        muF = lmbda["ineqnonlin"][:nmu]
        muT = lmbda["ineqnonlin"][nmu:nmu + nmu]
        if self.flow_lim == "I":
            dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It = \
                self.om.case.dIbr_dV(self._Yf, self._Yt, V)
            Hfaa, Hfav, Hfva, Hfvv = \
                self.om.case.d2AIbr_dV2(dIf_dVa, dIf_dVm, If, self._Yf, V, muF)
            Htaa, Htav, Htva, Htvv = \
                self.om.case.d2AIbr_dV2(dIt_dVa, dIt_dVm, It, self._Yt, V, muT)
        else:
            f = [e.from_bus._i for e in self._ln]
            t = [e.to_bus._i for e in self._ln]
            # Line-bus connection matrices.
            Cf = csr_matrix((ones(self._nl), (range(self._nl), f)), (self._nl, self._nb))
            Ct = csr_matrix((ones(self._nl), (range(self._nl), t)), (self._nl, self._nb))
            dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St = \
                self.om.case.dSbr_dV(self._Yf, self._Yt, V)
            if self.flow_lim == PFLOW:
                Hfaa, Hfav, Hfva, Hfvv = \
                    self.om.case.d2ASbr_dV2(dSf_dVa.real(), dSf_dVm.real(),
                                            Sf.real(), Cf, self._Yf, V, muF)
                Htaa, Htav, Htva, Htvv = \
                    self.om.case.d2ASbr_dV2(dSt_dVa.real(), dSt_dVm.real(),
                                            St.real(), Ct, self._Yt, V, muT)
            elif self.flow_lim == SFLOW:
                Hfaa, Hfav, Hfva, Hfvv = \
                    self.om.case.d2ASbr_dV2(
                        dSf_dVa, dSf_dVm, Sf, Cf, self._Yf, V, muF)
                Htaa, Htav, Htva, Htvv = \
                    self.om.case.d2ASbr_dV2(
                        dSt_dVa, dSt_dVm, St, Ct, self._Yt, V, muT)
            else:
                raise ValueError

        d2H = vstack([
            hstack([
                vstack([hstack([Hfaa, Hfav]),
                        hstack([Hfva, Hfvv])]) +
                vstack([hstack([Htaa, Htav]),
                        hstack([Htva, Htvv])]),
                csr_matrix((2 * self._nb, nxtra))
            ]),
            hstack([
                csr_matrix((nxtra, 2 * self._nb)),
                csr_matrix((nxtra, nxtra))
            ])
        ], "csr")

        return d2f + d2G + d2H