2015-07-27 14 views
1

Скажем, у меня есть линейная система с N переменных, но у меня есть только N-1 уравнений (ограничения). Как я могу получить выполнимый набор (диапазон) для каждой из N переменных, используя R?Оптимизация решения с линейными ограничениями

Пример:

A <- matrix(data=c(0,1,0,1,0,1,0,1, 
        0,0,1,1,0,0,1,1, 
        0,0,0,0,1,1,1,1, 
        0,0,0,1,0,0,0,1, 
        0,0,0,0,0,1,0,1, 
        0,0,0,0,0,0,1,1, 
        1,1,1,1,1,1,1,1 
        ), 
      ncol=8, byrow=T) 
b <- matrix(data=c(0.2,0.4,0.6,0.06,0.18,0.12,1), 
      ncol=1) 
> A 
##  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] 
##[1,] 0 1 0 1 0 1 0 1 
##[2,] 0 0 1 1 0 0 1 1 
##[3,] 0 0 0 0 1 1 1 1 
##[4,] 0 0 0 1 0 0 0 1 
##[5,] 0 0 0 0 0 1 0 1 
##[6,] 0 0 0 0 0 0 1 1 
##[7,] 1 1 1 1 1 1 1 1 
> b 
##  [,1] 
##[1,] 0.20 
##[2,] 0.40 
##[3,] 0.60 
##[4,] 0.06 
##[5,] 0.18 
##[6,] 0.12 
##[7,] 1.00 

Дополнительное ограничение: каждое значение должно быть possitive и должна быть в пределах от 0 до 1 (этот последний может быть опущен, так как последнее уравнение добавляет до 1 только с possitive значениями)

+0

Такого рода проблемы могут быть решены с помощью статистических инструментов регуляризован- регрессия. Один из них - знаменитый LASSO. –

ответ

1

Я бы рекомендовал использовать методы в качестве регуляризованной регрессии. Существуют пакеты в R, такие как lars или glmnet. Я бы посоветовал прочитать это answer.

В принципе, в R вы бы запустить

library(lars) 
model=lars(A,b) 
coef(model) 
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] 
[1,] 0.000 0.00 0.000 0 0.0000 0.0000 0.00000 0 
[2,] 0.000 0.00 0.000 0 0.0982 0.0000 0.00000 0 
[3,] 0.131 0.00 0.000 0 0.2000 0.0000 0.00000 0 
[4,] 0.138 0.00 0.185 0 0.3849 0.0000 0.00000 0 
[5,] 0.109 0.00 0.260 0 0.4000 0.0600 0.00000 0 
[6,] 0.109 0.00 0.260 0 0.3997 0.0603 0.00029 0 
[7,] 0.100 0.02 0.280 0 0.3600 0.1200 0.06000 0 

Результаты является матрица 7х8, где в каждой строке сверху вниз, наиболее важной переменной вводится в модель. Например, первая строка пытается подогнать модель только с перехватом. Вторая вводит самую сильную переменную 8, которая является 5-й, и так далее. Из результатов вы можете видеть, что переменные 4 и 8 все время имеют нулевой коэффициент, что означает, что они не важны, учитывая все остальные.

1

Вам необходимо вычислить 'nullspace' of A, например. с пакетом MASS или функция nullspace() in pracma. Во-первых вычисления по методу наименьших квадратов раствора, а затем добавить линейную комбинацию векторов в нуль-пространства:

library(pracma) 
N <- nullspace(A) 
# 0.3535534 * c(-1, 1, 1, -1, 1, -1, 1, -1) 
x0 <- qr.solve(A, b) 
# [1] 0.16 -0.04 0.22 0.06 0.30 0.18 0.12 0.00 

и х0 + х * Н, х реальный, генерирует все возможные решения.

N в этом примере является 1-мерным. Поскольку ранг A равен 1. При еще меньшем количестве уравнений пустое пространство будет иметь больше размеров.

+0

Есть ли способ добавить ограничение, что все решения должны быть приемлемыми? – Barranka

+0

Существует не такое x, что все компоненты x0 + x * N положительны, как вы можете проверить, выведя неравенства для x из всех компонентов (если я правильно вычислил). –

2

Проблема, которую вы описываете, вероятно, лучше всего подходит для линейного программирования, поскольку у вас есть набор линейных ограничений (Ax = b, x> = 0, x < = 1) и линейная цель (либо максимизировать, либо минимизировать некоторые x_i). Учитывая простую структуру вашей программы, я бы решить с помощью lpSolve пакета в R:

library(lpSolve) 
for (idx in seq_len(ncol(A))) { 
    for (type in c("min", "max")) { 
    mod <- lp(direction = type, 
       objective.in = as.numeric(seq_len(ncol(A)) == idx), 
       const.mat = rbind(A, diag(ncol(A))), 
       const.dir = c(rep("=", length(b)), rep("<=", ncol(A))), 
       const.rhs = c(b, rep(1, ncol(A)))) 
    print(paste("Variable:", idx, "type:", type, "value", mod$objval)) 
    } 
} 
# [1] "Variable: 1 type: min value 0.1" 
# [1] "Variable: 1 type: max value 0.12" 
# [1] "Variable: 2 type: min value 0" 
# [1] "Variable: 2 type: max value 0.02" 
# [1] "Variable: 3 type: min value 0.26" 
# [1] "Variable: 3 type: max value 0.28" 
# [1] "Variable: 4 type: min value 0" 
# [1] "Variable: 4 type: max value 0.02" 
# [1] "Variable: 5 type: min value 0.34" 
# [1] "Variable: 5 type: max value 0.36" 
# [1] "Variable: 6 type: min value 0.12" 
# [1] "Variable: 6 type: max value 0.14" 
# [1] "Variable: 7 type: min value 0.06" 
# [1] "Variable: 7 type: max value 0.08" 
# [1] "Variable: 8 type: min value 0.04" 
# [1] "Variable: 8 type: max value 0.06" 
Смежные вопросы