如何用R软件实现混合整数非线性规划(MINLP)的求解?另外有没有推荐的教材或典型的实例。谢谢!

Python013

如何用R软件实现混合整数非线性规划(MINLP)的求解?另外有没有推荐的教材或典型的实例。谢谢!,第1张

整数规划

integer programming

一类要求问题中的全部或一部分变量为整数的数学规划.

一般认为非线性的整数规划可分成线性部分和整数部分,因此常常把整数规划作为线性规划的特殊部分.在线性规划问题中,有些最优解可能是分数或小数,但对于某些具体问题,常要求解答必须是整数.例如,所求解是机器的台数,工作的人数或装货的车数等.为了满足整数的要求,初看起来似乎只要把已得的非整数解舍入化整就可以了.实际上化整后的数不见得是可行解和最优解,所以应该有特殊的方法来求解整数规划.在整数规划中,如果所有变量都限制为整数,则称为纯整数规划;如果仅一部分变量限制为整数,则称为混合整数规划.整数规划的一种特殊情形是01规划,它的变数仅限于0或1.

这里分享一下R语言实现VAR和SVAR的整个流程。

主要步骤包括:

1.单位根检验

2.确定滞后阶数

3.格兰杰因果检验

4.模型稳定性检验

5.脉冲响应

6.方差分解

(Johansen协整检验,如果需要的话)

整个过程用到的R语言的扩展包有:

library(zoo)

library(vars)

library(tseries)

首先,数据是下面的样子:

ps:数据是时间序列类型,可以通过下面方法将dataframe转成时间序列类型

data = ts(data)

1.单位根检验

#对data的第一列进行单位根检验

adf.test(data[,1])

2.滞后阶数确定

VARselect函数结果包括AIC、HQ、SC和FPE准则

#参数y为时间序列数据,lag.max为最大滞后阶数

#参数type值包括const截距,trend趋势,both同时包含截距和趋势,none不包含截距和趋势

VARselect(y=data, lag.max = 10, type = c("const"))

3.格兰杰因果检验

格兰杰因果检验有两个方法,第一个是在构造模型之前,第二个是在构造模型之后在模型的基础上进行格兰杰因果检验。

(1)构造模型之前格兰杰因果检验

#函数格式:grangertest(yt~xt)

eg:

grangertest(Value~BCI)

(2)构造模型之后格兰杰因果检验

#函数格式:causality(VARModel,cause)

eg

var =  VAR(data ,p = 2, type = "const")

causality(var,cause=c('Count','Value'))

ps:在这里如果想要构建SVAR模型的话,需要根据实际情况构建两个矩阵amat和bmat,然后使用这两个矩阵来构建SVAR模型:

svar = SVAR(var,Amat = amat,Bmat = bmat)

4.模型稳定性检验

#这里使用“OLS-CUSUM”,它给出的是残差累积和,在该检验生成的曲线图中,残差累积和曲线以时间为横坐标,

#图中绘出两条临界线,如果累积和超出了这两条临界线,则说明参数不具有稳定性。

sta = stability(var, type = c("OLS-CUSUM"), h = 0.15, dynamic = FALSE, rescale = TRUE)

plot(sta)##结果稳健

5.脉冲响应

#标题栏说明,这是BCI(或者其他变量)对各个变量(包括BCI自身)的脉冲响应

(1)VAR脉冲响应

var.irf<-irf(var,n.head=10)

plot(var.irf)

(2)SVAR脉冲响应

svar.irf<-irf(svar,n.ahead = 100)

plot(svar.irf)

6.方差分解

#反映了各变量的贡献率

(1)VAR方差分解

fevd1<-fevd(var, n.ahead = 10)

fevd1$Count

(2)SVAR方差分解

fevd2<-fevd(svar, n.ahead = 10)

fevd2$Value

ps:有时候需要进行Johansen协整检验

#Johansen协整检验,

#对r=0(不存在协整关系)的检验统计量大于临界值,表明拒绝原假设

yJoTest = ca.jo(data, type = c("trace"), ecdet = c("none"), K = 2)

summary(yJoTest)

网页链接

模型拟合

对于人口模型可以采用Logistic增长函数形式,它考虑了初期的指数增长以及总资源的限制。其函数形式如下。

首先载入car包以便读取数据,然后使用nls函数进行建模,其中theta1、theta2、theta3表示三个待估计参数,start设置了参数初始值,设定trace为真以显示迭代过程。nls函数默认采用Gauss-Newton方法寻找极值,迭代过程中第一列为RSS值,后面三列是各参数估计值。然后用summary返回回归结果。

library(car)

pop.mod1 <- nls(population ~ theta1/(1+exp(-(theta2+theta3*year))),start=list(theta1 = 400, theta2 = -49, theta3 = 0.025), data=USPop, trace=T)

summary(pop.mod)

在上面的回归过程中我们直接指定参数初始值,另一种方法是采用搜索策略,首先确定参数取值范围,然后利用nls2包的暴力方法来得到最优参数。但这种方法相当费时。

还有一种更为简便的方法就是采用内置自启动模型(self-starting Models),此时我们只需要指定函数形式,而不需要指定参数初始值。本例的logistic函数所对应的selfstarting函数名为SSlogis

pop.mod2 <- nls(population ~ SSlogis(year,phi1,phi2,phi3),data=USPop)

二、判断拟合效果

非线性回归模型建立后需要判断拟合效果,因为有时候参数最优化过程会捕捉到局部极值点而非全局极值点。最直观的方法是在原始数据点上绘制拟合曲线。

library(ggplot2)

p <- ggplot(USPop,aes(year, population))