Авторами построена осредненная численно-аналитическая теория движения второго порядка по массам планет для четырехпланетной задачи. Гамильтониан и уравнения движения даны в координатах Якоби и записаны в элементах второй системы Пуанкаре. Эксцентрические и облические элементы Пуанкаре сохранены в уравнениях движения вплоть третьей степени. В работе рассматривается орбитальная эволюция трехпланетной системы HD 39194 и четырехпланетных систем HD 141399, HD 160691 ( Ara). Численное интегрирование уравнений движения проводится для набора начальных условий, в которых неизвестные и известные из наблюдений с ошибками элементы орбит варьируются в допустимых пределах. Определены диапазоны изменения элементов орбит в зависимости от начальных условий. Предположение о стабильности наблюдаемых планетных систем позволяет исключить начальные условия, ведущие к экстремальному росту эксцентриситетов и наклонов орбит. Выявлены начальные условия, при которых элементы орбит сохраняются малыми на всем интервале моделирования. Показан способ, позволяющий сузить диапазон возможных значений неизвестных элементов орбит и определить их наиболее вероятные, с точки зрения стабильности, значения.