def solve(self):
        """ Runs dynamic simulation.

        @rtype: dict
        @return: Solution dictionary with the following keys:
                   - C{angles} - generator angles
                   - C{speeds} - generator speeds
                   - C{eq_tr} - q component of transient voltage behind
                     reactance
                   - C{ed_tr} - d component of transient voltage behind
                     reactance
                   - C{efd} - Excitation voltage
                   - C{pm} - mechanical power
                   - C{voltages} - bus voltages
                   - C{stepsize} - step size integration method
                   - C{errest} - estimation of integration error
                   - C{failed} - failed steps
                   - C{time} - time points
        """
        t0 = time()
        buses = self.dyn_case.buses

        solution = NewtonPF(self.case).solve()

        if not solution["converged"]:
            logger.error("Power flow did not converge. Exiting...")
            return {}
        elif self.verbose:
            logger.info("Power flow converged.")

        # Construct augmented Ybus.
        if self.verbose:
            logger.info("Constructing augmented admittance matrix...")

        gbus = [g.bus._i for g in self.dyn_generators]
        ng = len(gbus)

        Um = array([bus.v_magnitude for bus in buses])
        Ua = array([bus.v_angle * (pi / 180.0) for bus in buses])
        U0 = Um * exp(1j * Ua)
        U00 = U0

        augYbus = self.dyn_case.getAugYbus(U0, gbus)
        augYbus_solver = splu(augYbus)

        # Calculate initial machine state.
        if self.verbose:
            logger.info("Calculating initial state...")

        Efd0, Xgen0 = self.dyn_case.generatorInit(U0)
        omega0 = Xgen0[:, 1]

        Id0, Iq0, Pe0 = self.dyn_case.machineCurrents(Xgen0, U0)
        Vgen0 = r_[Id0, Iq0, Pe0]

        # Exciter initial conditions.
        Vexc0 = abs(U0[gbus])
        Xexc0, Pexc0 = self.dyn_case.exciterInit(Efd0, Vexc0)

        # Governor initial conditions.
        Pm0 = Pe0
        Xgov0, Pgov0 = self.dyn_case.governorInit(Pm0, omega0)
        Vgov0 = omega0

        # Check steady-state.
        Fexc0 = self.dyn_case.exciter(Xexc0, Pexc0, Vexc0)
        Fgov0 = self.dyn_case.governor(Xgov0, Pgov0, Vgov0)
        Fgen0 = self.dyn_case.generator(Xgen0, Xexc0, Xgov0, Vgen0)

        # Check Generator Steady-state
        if sum(abs(Fgen0)) > 1e-06:
            logger.error("Generator not in steady-state. Exiting...")
            return {}
        # Check Exciter Steady-state
        if sum(abs(Fexc0)) > 1e-06:
            logger.error("Exciter not in steady-state. Exiting...")
            return {}
        # Check Governor Steady-state
        if sum(abs(Fgov0)) > 1e-06:
            logger.error("Governor not in steady-state. Exiting...")
            return {}

        if self.verbose:
            logger.info("System in steady-state.")

        # Initialization of main stability loop.
        t = -0.02 # simulate 0.02s without applying events
        erst = False
        failed = False
        eulerfailed = False

        stoptime = self.dyn_case.stoptime

        if (isinstance(self.method, RungeKuttaFehlberg) or
            isinstance(self.method, RungeKuttaHighamHall)):
            stepsize = self.minstep
        else:
            stepsize = self.dyn_case.stepsize

        ev = 0
        eventhappened = False
        i = 0

        # Allocate memory for variables.
        if self.verbose:
            logger.info("Allocating memory...")

        chunk = 5000
        time = zeros(chunk)
        time[0, :] = t
        errest = zeros(chunk)
        errest[0, :] = erst
        stepsizes = zeros(chunk)
        stepsizes[0, :] = stepsize

        # System variables
        voltages = zeros(chunk)
        voltages[0, :] = U0.H

        # Generator
        angles = zeros((chunk, ng))
        angles[0, :] = Xgen0[:, 0] * 180.0 / pi
        speeds = zeros((chunk, ng))
        speeds[0, :] = Xgen0[:, 0] / 2 * pi * self.dyn_case.freq
        Eq_tr = zeros((chunk, ng))
        Eq_tr[0, :] = Xgen0[:, 2]
        Ed_tr = zeros((chunk, ng))
        Ed_tr[0, :] = Xgen0[:, 3]

        # Exciter and governor
        Efd = zeros((chunk, ng))
        Efd[0, :] = Efd0[:, 0]
        PM = zeros((chunk, ng))
        PM[0, :] = Pm0[:, 0]

        # Main stability loop.
        while t < stoptime + stepsize:
            i += 1
            if i % 45 == 0 and self.verbose:
                logger.info("%6.2f%% completed." % t / stoptime * 100)

            # Numerical Method.
            Xgen0, self.Pgen0, Vgen0, Xexc0, Pexc0, Vexc0, Xgov0, Pgov0, \
                Vgov0, U0, t, newstepsize = self.method.solve(t, Xgen0,
                    self.Pgen0, Vgen0, Xexc0, Pexc0, Vexc0, Xgov0, Pgov0,
                    Vgov0, augYbus_solver, gbus, stepsize)

#            if self.method == MODIFIED_EULER:
#                solver = ModifiedEuler(t, Xgen0, self.Pgen0, Vgen0, Xexc0,
#                                       Pexc0, Vexc0, Xgov0, Pgov0, Vgov0,
#                                       augYbus_solver, gbus, stepsize)
#
#                Xgen0, self.Pgen0, Vgen0, Xexc0, Pexc0, Vexc0, Xgov0, Pgov0,
#                Vgov0, U0, t, newstepsize = solver.solve()
#            elif self.method == RUNGE_KUTTA:
#                pass
#            elif self.method == RUNGE_KUTTA_FEHLBERG:
#                pass
#            elif self.method == HIGHAM_HALL:
#                pass
#            elif self.method == MODIFIED_EULER2:
#                pass
#            else:
#                raise ValueError

            if eulerfailed:
                logger.info("No solution found. Exiting... ")
                return {}

            if failed:
                t = t - stepsize

            # End exactly at stop time.
            if t + newstepsize > stoptime:
                newstepsize = stoptime - t
            elif stepsize < self.minstep:
                logger.info("No solution found with minimum step size. Exiting... ")
                return {}

            # Allocate new memory chunk if matrices are full.
            if i > time.shape[0]:
                time = zeros(chunk)
                errest = zeros(chunk)
                stepsize = zeros(chunk)
                voltages = zeros(chunk)
                angles = zeros((chunk, ng))
                speeds = zeros((chunk, ng))
                Eq_tr = zeros((chunk, ng))
                Ed_tr = zeros((chunk, ng))
                Efd = zeros((chunk, ng))
                PM = zeros((chunk, ng))

            # Save values.
            stepsizes[i, :] = stepsize
            errest[i, :] = erst
            time[i, :] = t

            voltages[i, :] = U0

            # Exciters
            Efd[i, :] = Xexc0[:, 0]
            # TODO: Set Efd to zero when using classical generator model.

            # Governors
            PM[i, :] = Xgov0[:, 0]

            # Generators
            angles[i, :] = Xgen0[:, 0] * 180.0 / pi
            speeds[i, :] = Xgen0[:, 1] * (2 * pi * self.dyn_case.freq)
            Eq_tr[i, :] = Xgen0[:, 2]
            Ed_tr[i, :] = Xgen0[:, 3]

            # Adapt step size if event will occur in next step.
            if (len(self.events) > 0 and ev <= len(self.events) and
                isinstance(self.method, RungeKuttaFehlberg) and
                isinstance(self.method, RungeKutta)):

                if t + newstepsize >= self.events[ev].t:
                    if self.events[ev] - t < newstepsize:
                        newstepsize = self.events[ev].t - t

            # Check for events.
            if len(self.events) > 0 and ev <= len(self.events):
                for event in self.events:
                    if (abs(t - self.events[ev].t) > 10 * EPS or
                        ev > len(self.events)):
                        break
                    else:
                        eventhappened = True

                    event.obj.set_attr(event.param, event.newval)

                    ev += 1

                if eventhappened:
                    # Refactorise.
                    self.dyn_case.getAugYbus(U00, gbus)
                    U0 = self.dyn_case.solveNetwork(Xgen0, self.Pgen0,
                                                    augYbus_solver, gbus)

                    Id0, Iq0, Pe0 = self.dyn_case.machineCurrents(Xgen0,
                                                                  self.Pgen0,
                                                                  U0[gbus])
                    Vgen0 = r_[Id0, Iq0, Pe0]
                    Vexc0 = abs(U0[gbus])

                    # Decrease stepsize after event occured.
                    if (isinstance(self.method, RungeKuttaFehlberg) or
                        isinstance(self.method, RungeKuttaHighamHall)):

                        newstepsize = self.minstepsize

                    # If event occurs, save values at t- and t+.
                    i += 1

                    # Save values
                    stepsize[i, :] = stepsize.T
                    errest[i, :] = erst.T
                    time[i, :] = t

                    voltages[i, :] = U0.T

                    # Exciters.
                    # Set Efd to zero when using classical generator model.
#                    Efd[i, :] = Xexc0[:, 1] * (flatnonzero(genmodel > 1))

                    # Governors.
                    PM[i, :] = Xgov0[:, 1]

                    # Generators.
                    angles[i, :] = Xgen0[:, 0] * 180.0 / pi
                    speeds[i, :] = Xgen0[:, 1] / (2.0 * pi * self.freq)
                    Eq_tr[i, :] = Xgen0[:, 2]
                    Ed_tr[i, :] = Xgen0[:, 3]

                    eventhappened = False

            # Advance time
            stepsize = newstepsize
            t += stepsize

        # End of main stability loop ------------------------------------------

        # Output --------------------------------------------------------------

        if self.verbose:
            logger.info("100%% completed")
            elapsed = time() - t0
            logger.info("Simulation completed in %5.2f seconds." % elapsed)

        # Save only the first i elements.
        angles = angles[0:i, :]
        speeds = speeds[0:i, :]
        Eq_tr = Eq_tr[0:i, :]
        Ed_tr = Ed_tr[0:i, :]

        Efd = Efd[0:i, :]
        PM = PM[0:i, :]

        voltages = voltages[0:i, :]

        stepsize = stepsize[0:i, :]
        errest = errest[0:i, :]
        time = time[0:i, :]

        if self.plot:
            raise NotImplementedError

        return {}