We consider region of attraction for a dynamical system - a set of points from which system trajectories go to a stable equilibrium. One can estimate it via the sublevel sets of a positive definite Lyapunov function, which is decreasing along any trajectory in the region. We formulate this task as a bilevel optimization over the function parameters and system state variables. On the lower level we have minimization of the function with nonconvex equality constraint. On the upper level a size estimation of a sublevel set needs to be maximized. Local improvement algorithm is discussed.