Я пытаюсь решить систему из n уравнений в R. n-1 уравнений нелинейны, а последнее линейно. Если это поможет, это задача оптимизации с ограничениями, где уравнение n-1 является условием первого порядка, а последнее — бюджетным ограничением. Нелинейные уравнения n-1 можно охарактеризовать: нелинейными уравнениями
Если изображение не отображается, его можно определить поэлементно, например:
v_i*epsilon_i*cos(2/pi * e_i/delta_i)-лямбда=0
где эпсилон, v, e и дельта — векторы размерности n-1, а лямбда — скаляр, общий для всех уравнений).
Последнее уравнение представляет собой простое линейное уравнение вида |e|=c. То есть сумма всех элементов в e есть некоторое известное c, называемое ниже parms[1,4] или «бюджет».
Меня интересует решение для вектора e и постоянной лямбды, рассматривая все остальное как данность.
Я попытался использовать мультирут rootSolve. Для этого я определяю единственный вектор X, который должен быть вектором e, с добавленной к нему лямбдой, так что многокорень решает для 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,])
Когда я запускаю это, я получаю сообщение об ошибке
Ошибка в stode(y, times, func, parms = parms, ...): REAL() может применяться только к «числовому», а не к «списку»
Чего я действительно не понимаю - как список функций может не относиться к классу "список"? Итак, мой второй вопрос:
- Как создать список функций, если они не являются простыми однострочными функциями (как в ссылках выше)
Наконец, я был бы очень признателен за любые рекомендации о том, как лучше решать эти типы проблем в R.
Спасибо за любую помощь!
