- 积分
- 7
- 注册时间
- 2002-9-10
- 仿真币
-
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2009-12-6 20:50:47
|
显示全部楼层
来自 山东淄博
本帖最后由 wanglu 于 2009-12-7 17:32 编辑
下面我把我的问题发出来，希望大家能帮忙了。。。
napler821127 发表于 2009-12-4 16:05 
napler821127的问题，Forcal的优化效果也不好。经napler821127验证，Matlab的ODE求解器与Forcal使用的“积分一步的连分式法”对微分方程的计算结果一致，故优化效果不理想可能来自“求n维极值的单形调优法”，但目前Forcal中还没有更好的优化方法。
故我想到使用Matlab与FcScript（Forcal脚本服务器）混合编程求解：使用Matlab的优化函数进行参数优化，微分方程计算部分由Forcal完成，优化结果由Matlab验证。虽然接口部分效率较低，但Forcal的快速计算速度足以抵消这个不足，故总的解题时间预计可以减少。
FcScript是由Forcal和MForcal支持的脚本控件，采用双接口、STA单线程套间模型。目前定义了两种接口：VBMForcal接口（用于VBScript脚本，全部使用VARIANT参数）和CMForcal接口（用于Excel、CAD等各种使用COM对象的地方）。
FcScript下载：http://xoomer.virgilio.it/forcal/xiazai/forcal9/fcscript.rar
napler821127的问题描述：
为进行混合编程，针对本问题，以求解300000时的数据为例，Forcal提供了两组函数，每组函数编译为一个模块，各提供一个函数供Matlab调用。Matlab根据实际情况选择其中的一个即可。模块1中提供了函数optint，该函数根据优化参数用连分式法对微分方程组进行积分；模块2中提供了函数opt，直接根据优化参数返回目标函数值。
模块1、用连分式法对微分方程组进行积分，获得理论值（可简称为optint）
-
- !using["XSLSF"]; //使用命名空间XSLSF
- //函数定义，用于计算微分方程组中各方程右端函数值，连分式法对微分方程组积分一步函数pbs1将调用该函数f。
- //t为自变量，f1,f2为函数值，df1,df2为右端函数值（即微分方程的值）。
- //该函数被调用时将用到参数a1,a2,a3,a4,a5，这5个参数正好是拟合参数，用模块变量传递这5个参数。
- f(t,f1,f2,df1,df2 : x : a1,a2,a3,a4,a5)=
- {
- x=300000,
- df1=a1*(x-f2)^a2,
- df2=a3*df1-a4*f2^a5
- };
- //用连分式法对微分方程组进行积分，获得理论值。用编译符“~”输出该函数，供Matlab等客户机程序调用。
- //_a1,_a2,_a3,_a4,_a5为优化参数；t1,t2为积分的起点和终点；f1,f2为t1时的初值，返回时存放t2时的积分值。
- //函数返回时，t1=t2，或者t1-t2是一个极小值。下一次积分时，用t1,f1,f2作为积分起点比用t2,f1,f2要好。
- //h,s为自动变量。Array,free,hf,step,eps为静态变量：Array为工作数组；hf为函数f的句柄；step为积分步长；eps为积分精度。
- //专用静态变量free用于销毁表达式前的释放工作（Forcal在销毁表达式前将自动设置free=1，然后自动执行表达式）。
- ~optint(_a1,_a2,_a3,_a4,_a5,t1,f1,f2,t2 : h,s,static,Array,free,hf,step,eps : a1,a2,a3,a4,a5)=
- {
- if[free,delete(Array),return(0)], //销毁表达式时，销毁工作数组Array
- if{ !hf,
- hf=HFor("f"), //预先用函数HFor获得函数f的句柄hf
- Array=new[rtoi(real_s),rtoi(2*15)], //预先申请工作数组Array
- step=0.1,eps=1e-6 //积分步长step要合适，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1
- },
- a1=_a1,a2=_a2,a3=_a3,a4=_a4,a5=_a5,//传递优化变量，函数f中要用到a1,a2,a3,a4,a5
- Array.setra(0,f1,f2), //设置积分初值，通过模块变量Array传递，Array是一个数组
- s=ceil[(t2-t1)/step], //计算积分步数
- h=(t2-t1)/s, //重新计算积分步长
- { pbs1[hf,t1,Array,h,eps], //对微分方程组积分一步
- t1=t1+h //积分步长增加
- }.until[abs(t1-t2)<h/2], //连续积分至t2
- Array.getra(0,&f1,&f2) //从数组Array获得t2时的积分值
- };
- !HFor("optint"); //使函数optint有效
模块2、获得优化时目标函数的值（可简称为opt）
-
- !using["XSLSF"]; //使用命名空间XSLSF
- //函数定义，用于计算微分方程组中各方程右端函数值，连分式法对微分方程组积分一步函数pbs1将调用该函数ff。
- //t为自变量，f1,f2为函数值，df1,df2为右端函数值（即微分方程的值）。
- //该函数被调用时将用到参数a1,a2,a3,a4,a5，这5个参数正好是拟合参数，用模块变量传递这5个参数。
- ff(t,f1,f2,df1,df2 : x : a1,a2,a3,a4,a5)=
- {
- x=300000,
- df1=a1*(x-f2)^a2,
- df2=a3*df1-a4*f2^a5
- };
- //用连分式法对微分方程组进行积分，获得理论值。
- //t1,t2为积分的起点和终点。
- //h,s为自动变量。
- //模块变量：hff为函数ff的句柄，要预先获得该句柄；Array为工作数组；step为积分步长；eps为积分精度。
- 积分(t1,t2:h,s:hff,Array,step,eps)=
- {
- s=ceil[(t2-t1)/step], //计算积分步数
- h=(t2-t1)/s, //重新计算积分步长
- { pbs1[hff,t1,Array,h,eps], //对微分方程组积分一步
- t1=t1+h //积分步长增加
- }.until[abs(t1-t2)<h/2] //连续积分至t2
- };
- //获得优化时目标函数的值。用编译符“~”输出该函数，供Matlab等客户机程序调用。
- //_a1,_a2,_a3,_a4,_a5为优化参数；t1,t2,i,s,y,yy,z为自动变量；tArray,free,max为静态变量；hff,Array,step,eps, a1,a2,a3,a4,a5为模块变量。
- //数组tArray中存放实验数据；Array为工作数组；hff为函数ff的句柄；step为积分步长；eps为积分精度。
- //专用静态变量free用于销毁表达式前的释放工作（Forcal在销毁表达式前将自动设置free=1，然后自动执行表达式）。
- ~opt(_a1,_a2,_a3,_a4,_a5 : t1,t2,i,s,y,yy,z, static,tArray,free,max : hff,Array,step,eps, a1,a2,a3,a4,a5)=
- {
- if[free,delete(Array),delete(tArray),return(0)], //销毁表达式时，销毁数组Array和tArray
- if{ !tArray,
- hff=HFor("ff"), //预先用函数HFor获得函数ff的句柄hff
- Array=new[rtoi(real_s),rtoi(2*15)], //预先申请工作数组Array
- max=4801, //实验数据组数
- tArray=new{rtoi(real_s),rtoi(max),rtoi(2),rtoi(EndType), //存放实验数据ti
- 0 , 0 ,
- 0.5 , 0.000809 ,
- 1 , 0.001182 ,
- //其余数据略
- 2399.5 , 0.006192 ,
- 2400 , 0.006192
- },
- step=0.1,eps=1e-6 //积分步长step要合适，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1
- },
- a1=_a1,a2=_a2,a3=_a3,a4=_a4,a5=_a5,//传递优化变量，函数ff中要用到a1,a2,a3,a4,a5
- Array.setra(0,0,50000), //设置积分初值，通过模块变量Array传递，Array是一个数组
- i=1,s=0,t1=0,
- (i<max).while{
- tArray.getra(i*2,&t2,&yy),
- 积分(&t1,t2), //从t1积分到t2
- Array.getra(0,&y,z), //从数组Array获得t2时的积分值
- s=s+[y-yy]^2, //累加理论值与实验值差的平方
- i++
- },
- s
- };
- !HFor("opt"); //使函数opt有效
现在的问题是如何在Matlab中调用函数optint或opt？望各位高手帮忙。
以下给出VBScript使用VBMFocal接口的例子及VBA使用CMForcal的例子，供参考。
1、VBScript使用VBMFocal的例子（步骤见注释部分）：调用模块1中的函数optint
-
- Set obj=CreateObject("FcScript.VBMForcal")
- dim str_optint, nModule, hModule, err1, err2, iErrCode, hHandle, NumPara, xx
- Call obj.CLoadDll("FcData32W","XSLSF32W") ‘加载需要的Forcal扩展动态库
- nModule = 0 ‘设为0
- str_optint=“模块2的代码” ‘用实际代码取代之
- iErrCode=obj.VBComModule (str_optint, nModule, hModule, err1, err2 ) '编译Forcal源程序为一个模块
- If iErrCode=0 Then
- hHandle=obj.VBGetForHandle(CLng(2)," optint ",NumPara)'获得函数optint的句柄
- xx =Array(_a1,_a2,_a3,_a4,_a5,t1,f1,f2,t2) '用数组xx存放自变量，设为有意义的实际值
- obj.VBCalFor(CLng(2), hHandle, xx) '执行函数optint，结果在xx中，可多次调用该函数
- Else
- MsgBox iErrCode,0,"编译错误代码"
- End If
- Set obj=Nothing
-
- Dim obj As Object
- Dim iCode As Long, theStr As String, nModule As Long, hModule As Long, err1 As Long, err2 As Long
- Set obj = CreateObject("FcScript.CMForcal")‘获得接口
- Call obj.CLoadDll("FcData32W","XSLSF32W") ‘加载需要的Forcal扩展动态库
- nModule = 0 ‘设为0
- theStr = “模块2的代码” ‘用实际代码取代之
- iCode = obj.CComModule(theStr, nModule, hModule, err1, err2) ‘编译模块2
- If iCode = 0 Then ‘检查编译是否正确
- ExComModule = "编译正确！" & "模块句柄：" & hModule
- Else
- ExComModule = "编译错误！错误代码：" & iCode & "，出错起始位置：" & err1 & "，出错结束位置：" & err2
- End If
- Dim hFor As Long, Para As Long, NumPara As Long
- hFor = obj.CGetForHandle(2, “opt”, Para, NumPara) ‘获得函数opt的句柄
- Dim d(5) As Double
- d(0) = 0.0 ‘设置优化参数
- d(1) = 1.0
- d(2) = 2.0
- d(3) = 3.0
- d(4) = 4.0
- Dim ExCalRFor As Double
- ExCalRFor = obj.CCalRFor(hFor, d(0)) ‘传入优化参数，获得目标函数值，可多次调用该函数
- Call obj.CDeleteModule(hModule) ‘销毁模块2
- Set obj = Nothing
以上虽然是用Matlab为例进行说明的，但有优化函数且能使用COM组件的软件均可一试。欢迎大家批评指正。 |
本帖子中包含更多资源
您需要 登录 才可以下载或查看，没有账号？注册
×
评分
-
1
查看全部评分
-
|
IP卡
狗仔卡
发表于 2009-12-4 16:05:32
显身卡
发表于 2009-12-7 21:53:11