Как решить n нелинейных уравнений в R на основе векторов / матриц - PullRequest
0 голосов
/ 04 мая 2020

Я пытаюсь решить систему из n уравнений в R. n-1 уравнений нелинейны, а последнее - линейно. Если это помогает, это ограниченная задача оптимизации, где уравнение n-1 является условиями первого порядка, а последнее - бюджетным ограничением. Н-1 нелинейные уравнения могут быть охарактеризованы следующим образом: нелинейные уравнения

Если изображение не отображается, его можно определить поэлементно, например:

v_i * epsilon_i * cos (2 / pi * e_i / delta_i) -lambda = 0

, где epsilon, v, e и delta являются векторами n- 1 измерение, и лямбда является скаляром, общим для всех уравнений).

Последнее уравнение является простым линейным уравнением вида | e | = c. То есть сумма всех элементов в e - это некоторая известная c, называемая ниже parms [1,4] или "бюджет".

Я заинтересован в решении для вектора e и постоянной лямбды, относиться ко всему остальному как данному.

Я пытался использовать rootSolve's multi root. Чтобы сделать это, я определяю один вектор X, который должен быть вектором e, с добавлением лямбды к нему, так что multi root решает для x, и получает n уравнений в виде списка. Все параметры сохраняются в матрице под названием «parms»

Сначала я определяю n-1 нелинейных уравнений

convex_focs <- function(x = numeric(),parms = numeric()){
deltas = parms[,1]
v = parms[,2]
lambda = x[1]
e = x[2:length(x)]
epsilon_2 = exp(parms[,3]) - parms[,1]
return(epsilon_2*cos((pi/2)*(e/deltas))-lambda)
}

В этом уравнении используется матричная запись, и она отлично работает сама по себе , Затем я определяю последнее линейное уравнение:

convex_budget <- function(x = numeric(),parms = numeric()){
e = x[2:length(x)]
return(sum(e)-parms[1,4])
}

Затем я пробовал convex_system <- function(x,parms) c(convex_focs,convex_budget ) и вызываю:

multiroot(f = convex_system, maxiter = 500000, start =  c(0,rep(budget/length(parms[,1]),length(parms[,1]))), parms = parms[,])

Это, конечно, не работает, так как rootSolve распознает convex_system как два уравнения, но X как n-мерность.

Если я опущу последнее уравнение и трактую лямбду как заданную (таким образом, только решая нелинейные уравнения) , я могу получить решение . Но, конечно, это нехорошо, потому что я не знаю лямбду.

Итак, мой первый вопрос: 1. Как мне создать список функций из моих векторов, которые rootSolve распознает как систему? Я пытался использовать lapply или виньетки, чтобы создать список convex_focs уравнений, по одному на каждый элемент в векторе, но не мог придумать, как заставить это работать. 2. Почему он распознал мою исходную функцию convex_focs как систему уравнений, но когда я добавил convex_budget, она перестала работать?

Я тогда (в отчаянии ...) попытался вручную определить набор функции, глядя только на 3 нелинейных, а не n-1. Это связано с тем, что мой список функций будет выглядеть как руководство и другие решения, которые я нашел в Интернете :

convex_system <- function(x,parms) c(F1 = 
                                     function(x =x,parms = parms){
                                       deltas = parms[1,1]
                                       v = parms[1,2]
                                       lambda = x[1]
                                       e = x[2]
                                       epsilon_2 = exp(parms[1,3]) - parms[1,1]
                                       return(v*epsilon_2*cos((pi/2)*(e/deltas))-lambda)
                                     }
                                     ,
                                   F2 = 
                                     function(x = x,parms = parms){
                                       deltas = parms[2,1]
                                       v = parms[2,2]
                                       lambda = x[1]
                                       e = x[3]
                                       epsilon_2 = exp(parms[2,3]) - parms[2,1]
                                       return(v*epsilon_2*cos((pi/2)*(e/deltas))-lambda)
                                     }
                                   ,
                                   F3 = 
                                     function(x = x,parms = parms){
                                       deltas = parms[3,1]
                                       v = parms[3,2]
                                       lambda = x[1]
                                       e = x[4]
                                       epsilon_2 = exp(parms[3,3]) - parms[3,1]
                                       return(v*epsilon_2*cos((pi/2)*(e/deltas))-lambda)
                                     }
                                     ,
                                   F_budget = function(x = x,parms = parms){
                                       e = x[2:length(x)]
                                       return(sum(e)-parms[1,4])}
                                  )

И вызов multiroot(f = convex_system, maxiter = 500000, start = c(0,rep(budget/length(parms[1:3,1]),length(parms[1:3,1]))), parms = parms[1:3,]) Когда я запускаю это Я получаю сообщение об ошибке

Ошибка в краже (y, times, fun c, parms = parms, ...): REAL () может быть применена только к 'цифре c' а не «список»

Что я действительно не понимаю - как может список функций не относиться к классу «список»? Итак, мой второй вопрос:

Как создать список функций, если они не являются простыми однострочными функциями (как в приведенных выше ссылках)

Наконец, я очень признателен за любые указания о том, как Лучше решить такие проблемы в R.

Спасибо за любую помощь!

1 Ответ

0 голосов
/ 05 мая 2020

Основной проблемой вашего подхода является спецификация функции convex_system. То, как вы это написали, подразумевает, что это вектор функций, и он не оценивается. Просто попробуйте один оператор convex_system(start,parms), чтобы увидеть возвращаемое значение.

Измените его на

convex_system <- function(x,parms) c(convex_focs(x,parms),convex_budget(x,parms) )

, который возвращает значения, возвращенные для указанных c значений x и parms.

Вы не предоставили никаких значений для констант и переменных. Поэтому мы не можем что-то попробовать.

Так что используйте поддельные данные:

budget <- 100
lambda <- 5

parms <- matrix(c(1,2,3,
                  2,3,4,
                  3,4,5,
                  105,0,0), ncol=4,byrow=FALSE)

parms
xstart <- c(0,rep(budget/length(parms[1:3,1]),length(parms[1:3,1])))

И, пожалуйста, не забудьте показать весь соответствующий код, даже library операторов.

Я попробовал два пакета для решения системы нелинейных уравнений: nleqslv и rootSolve.

library(nleqslv)
nleqslv(xstart,convex_system,parm=parms)

, что привело к

$x
[1] -18.07036  37.79143  34.44652  32.76205

$fvec
[1]  6.578382e-10 -4.952128e-11 -1.673328e-12  0.000000e+00

$termcd
[1] 1

$message
[1] "Function criterion near zero"

$scalex
[1] 1 1 1 1

$nfcnt
[1] 146

$njcnt
[1] 7

$iter
[1] 92

См. Документацию nleqslv для смысл элементов вышеприведенного списка. Кроме того, nleqslv в этом случае использовал метод Бройдена, который сохраняет вычисления Якобиана.

Использование rootSolve дает:

library(rootSolve)
multiroot(f = convex_system, maxiter = 500000, start = xstart, parms=parms)
$root
[1] -18.07036  37.79143  34.44652  32.76205

$f.root
[1] 1.650808e-06 4.383290e-08 8.365979e-08 1.250555e-11

$iter
[1] 10

$estim.precis
[1] 4.445782e-07

Как вы можете видеть, оба дают одинаковые результаты, но результаты nleqslv кажется ближе к нулю для значений составляющих функций (сравните fvec с f.root). Вы должны знать о разнице в критериях конвергенции (см. Документацию).

Решит ли это всю вашу проблему, зависит от вас.

Приложение

Кажется, что nleqslv нужно больше итераций, чем rootSolve. Это связано с используемым методом глобального поиска. Используя функцию testnslv, можно искать глобальный метод с использованием меньшего числа итераций, подобных этому

testnslv(xstart,convex_system,parm=parms)

со следующим результатом

Call:
testnslv(x = xstart, fn = convex_system, parm = parms)

Results:
    Method Global termcd Fcnt Jcnt Iter Message     Fnorm
1   Newton  cline      1   21   12   12   Fcrit 1.770e-24
2   Newton  qline      1   21   12   12   Fcrit 2.652e-24
3   Newton  gline      1   17    8    8   Fcrit 5.717e-25
4   Newton pwldog      1   45   31   31   Fcrit 9.837e-24
5   Newton dbldog      1   34   26   26   Fcrit 1.508e-25
6   Newton   hook      1   65   40   40   Fcrit 2.025e-26
7   Newton   none      1   10   10   10   Fcrit 7.208e-25
8  Broyden  cline      1   19    2   12   Fcrit 1.775e-19
9  Broyden  qline      1   19    2   12   Fcrit 1.768e-19
10 Broyden  gline      1   43    3   13   Fcrit 9.725e-18
11 Broyden pwldog      1  161    4  105   Fcrit 1.028e-19
12 Broyden dbldog      1  168    5  111   Fcrit 9.817e-21
13 Broyden   hook      1  121    7   67   Fcrit 5.138e-25
14 Broyden   none      1   11    1   11   Fcrit 7.487e-22

Можно видеть, что для method="Newton" методы gline и none (чистый Ньютон-Рафсон) требуют наименьшего количества итераций. И что метод Бройдена без глобального поиска вообще требует наименьшего количества оценок функций.

Предупреждение

Чтобы понять, почему некоторые из глобальных методов "лучше" укажите control=list(trace=1) в качестве аргумента nleqslv, например, global="none" и global="gline". Вы увидите, что чистый Ньютон не уменьшает критерий функции на каждой итерации. Просто повезло.

...