Я решаю задачу оптимизации с помощью пакета nsga2 в R, чтобы получить фронт Парето.
Я хочу минимизировать функцию, которая имеет в качестве входного вектора вектор из двух компонентов и имеет выходное значение от 0 до 3 ок.В то же время мне нужно максимизировать функцию с тем же входным значением и значением выходного сигнала от 0 до 6. Это код:
library("nsga2R")
library("deSolve")
library("ggplot2")
library("tidyverse")
library("bioinactivation")
library("dplyr")
library("tidyr")
value_inactivation <- function(x) {
time1 <- x[1]
temp <- x[2]
model_parms <- c(delta_ref = 3.9,
z = 4.2,
p = 1.2,
temp_ref = 55,
N0=10^6
)
times <- seq(0,time1, by = 10^-1)
temp_profile <- data.frame(time = c(0, time1), temperature = c(temp, temp))
prediction <- predict_inactivation("Mafart", times, model_parms, temp_profile)
a <- prediction$simulation %>%
filter(time == time1) %>%
mutate(res = logS) %>%
select(res) %>%
mutate(res = ifelse(res == "NaN" || res < -6, 6, -res)) %>%
mutate(res = res/6)
a <- a$res
return(a)
}
value_Maillard <- function(x) {
time1 <- x[1]
tempe <- x[2]
comp <- "Acrylamide"
a <- data.frame(temp = c(120,140,160,180,200),
k1 = c(0.55*10^-3,1.4*10^-3,3.1*10^-3,6.5*10^-3,13*10^-3),
k2 = c(0.02*10^-3, 0.078*10^-3, 0.27*10^-3, 0.86*10^-3, 2.4*10^-3),
k3 = c(1.9*10^-3 ,3.1*10^-3,5*10^-3,7.8*10^-3,12*10^-3),
k5 = c(30*10^-3,55*10^-3,95*10^-3,157*10^-3,246*10^-3),
k7 = c(1.7*10^-3,5.2*10^-3,14*10^-3,36*10^-3,84*10^-3),
k8 = c(3.9*10^-3,8.3*10^-3,16*10^-3,30*10^-3,53*10^-3),
k9 = c(2*10^-3,7.1*10^-3,22*10^-3,62*10^-3,159*10^-3))
#ggplot(a) + geom_line(aes(temp,log(k1))) +geom_line(aes(temp,log(k2))) +
# geom_line(aes(temp,log(k3))) +
# geom_line(aes(temp,log(k5))) +
# geom_line(aes(temp, log(k7))) +
# geom_line(aes(temp, log(k8))) +
# geom_line(aes(temp, log(k9)))
##Funcion log(k) = B0_k * temp + B1_k
asa <- summary(lm(data = a,formula = log(k1) ~ temp))$coefficients
B1_k1 <- asa[[1]]
B0_k1 <- asa[[2]]
R2_k1 <- summary(lm(data = a,formula = log(k1) ~ temp))$r.squared
asa <- summary(lm(data = a,formula = log(k2) ~ temp))$coefficients
B1_k2 <- asa[[1]]
B0_k2 <- asa[[2]]
R2_k2 <- summary(lm(data = a,formula = log(k2) ~ temp))$r.squared
asa <- summary(lm(data = a,formula = log(k3) ~ temp))$coefficients
B1_k3 <- asa[[1]]
B0_k3 <- asa[[2]]
R2_k3 <- summary(lm(data = a,formula = log(k3) ~ temp))$r.squared
asa <- summary(lm(data = a,formula = log(k5) ~ temp))$coefficients
B1_k5 <- asa[[1]]
B0_k5 <- asa[[2]]
R2_k5 <- summary(lm(data = a,formula = log(k5) ~ temp))$r.squared
asa <- summary(lm(data = a,formula = log(k7) ~ temp))$coefficients
B1_k7 <- asa[[1]]
B0_k7 <- asa[[2]]
R2_k7 <- summary(lm(data = a,formula = log(k7) ~ temp))$r.squared
asa <- summary(lm(data = a,formula = log(k8) ~ temp))$coefficients
B1_k8 <- asa[[1]]
B0_k8 <- asa[[2]]
R2_k8 <- summary(lm(data = a,formula = log(k8) ~ temp))$r.squared
asa <- summary(lm(data = a,formula = log(k9) ~ temp))$coefficients
B1_k9 <- asa[[1]]
B0_k9 <- asa[[2]]
R2_k9 <- summary(lm(data = a,formula = log(k9) ~ temp))$r.squared
# if(R2_k1 > 0.99 && R2_k2 > 0.99 && R2_k3 > 0.99 && R2_k5 > 0.99
# && R2_k7 > 0.99 && R2_k8 > 0.99) {
# print("El ajuste es bueno")
# } else { print("Mal ajuste")}
## Parameters k values
k1 <- function(tempe) {
k1 = exp(B0_k1 * tempe + B1_k1)
return(k1)
}
k2 <- function(tempe) {
k2 = exp(B0_k2 * tempe + B1_k2)
return(k2)
}
k3 <- function(tempe) {
k3 = exp(B0_k3 * tempe + B1_k3)
return(k3)
}
k5 <- function(tempe) {
k5 = exp(B0_k5 * tempe + B1_k5)
return(k5)
}
k7 <- function(tempe) {
k7 = exp(B0_k7 * tempe + B1_k7)
return(k7)
}
k8 <- function(tempe) {
k8 = exp(B0_k8 * tempe + B1_k8)
return(k8)
}
k9 <- function(tempe) {
k9 = exp(B0_k9 * tempe + B1_k9)
return(k9)
}
closed.sir.model <- function (t, x, tempe) {
##first extract the state variables
X1 <- x[1]
X2 <- x[2]
X3 <- x[3]
X4 <- x[4]
X5 <- x[5]
X6 <- x[6]
X7 <- x[7]
X8 <- x[8]
X9 <- x[9]
##Extract the parameters
k1 <- k1(tempe)
k2 <- k2(tempe)
k3 <- k3(tempe)
k5 <- k5(tempe)
k7 <- k7(tempe)
k8 <- k8(tempe)
k9 <- k9(tempe)
#Code the model equations
dX1 <- k9*X3
dX2 <- k8*X5
dX3 <- k1*X4-k9*X3
dX4 <- -k1*X4-k7*X4-k2*X4*X5
dX5 <- -k2*X4*X5-k8*X5
dX6 <- k7*X4
dX7 <- k2*X4*X5-k5*X7-k3*X7
dX8 <- k3*X7
dX9 <- k3*X7
##combine results into a single vector
dx <- c(dX1,dX2,dX3,dX4,dX5,dX6,dX7,dX8,dX9)
##return result as a list
list(dx)
}
##params <-c(k1=13*10^-3,k2=2.4*10^-3,k3=12*10^-3,k5=246*10^-3,k7=84*10^-3,k8=53*10^-3,k9=159*10^-3)
#params <- ks[[tempe]]
times <- seq(from=0, to=time1, by=1e-3)
xstart <- c(X1=0, X2=0, X3=0, X4=100, X5=100, X6=0, X7=0, X8=0, X9=0)
out <- as.data.frame(
ode(
func=closed.sir.model,
y=xstart,
times = times,
parms = tempe
)
)
out <- out %>% filter(time == time1)
switch(comp,
Formic_acid={
b <-out[["X1"]]
},
X2={
b <-out[["X2"]]
},
Glucose={
b <-out[["X3"]]
},
Fructose={
b <-out[["X4"]]
},
Aspargine={
b <-out[["X5"]]
},
Acetic_acid={
b <-out[["X6"]]
},
Schiff_base={
b <-out[["X7"]]
},
Acrylamide={
b <-out[["X8"]]
},
Melanoidins={
b <-out[["X9"]]
},
)
return(b)
}
f <- function(x) {value_Maillard(x)}
g <- function(x) {value_inactivation(x)}
nsga2(f,idim = 2, odim= 1,lower.bounds = c(0.1,50), upper.bounds = c(20, 100), constraints = g, cdim=1)
Я ожидаю фронт Парето в качестве вывода, но фактическим выводом является следующая ошибка:
Оценка векторизованной целевой функции вернула вектор.Вместо этого вы должны вернуть матрицу 1 на 100!