AMPL Matlab 自动求导（AD）求解雅克比、海森矩阵

1.Motivation

最非线性问题研究，不可避免的需要求解雅克比，海森矩阵，这部分工作繁琐，易出错，费时。幸运的AMPL不仅具有快速建模的功能，同时还可以给出雅克比和海森矩阵，最近花了一些时间对AMPL建立出来的模型在matlab下求解雅克比和海森矩阵信息的软件库进行编译，相信这项工作可以大幅度减少以后研究工作的工作量。

2. AMPL简介

AMPL（A Mathematical Programming Language）是一种描述并求解大规模复杂数学问题的建模语言。AMPL支持世界上大部分的求解器，如CBC、CPLEX、FortMP、Gurobi、MINOS、IPOPT、SNOPT和KNITRO。AMPL语言的的一个主要的特点是对优化问题的数学表达式的简化，这使得简明地可读地定义优化问题成为可能。根据NEOS的统计AMPL是使用最为广泛的数学模型语言。
同时AMPL提供了用户自定义求解器的接口，用户可以对建立好的模型开发自己的求解器，这是AMPL最厉害的地方，这也是相对于GAMS，我选择AMPL的原因。

3. 准备资料

(我这列出我成功编译所需的软件)
Visual studio
Matlab 2012b(32位)
AMPL 在github 上的开源包：http://ampl.com/resources/hooking-your-solver-to-ampl/ https://github.com/ampl/mp

4. 编译

4.1 编译amplsolv.lib文件

将http://ampl.com/netlib/ampl/solvers/index.html 全部下载，该文件夹下包含makefile 文件，其中makefile.vc是针对visual studio编写的。打开Visual studio命令行，cd 到该目录下，输入nmake /f makefile.vc 会在目录下编译出amplsolv.lib。在C++/C语言，头文件加这个lib文件足够了。但是我们是在matlab下运行，还需要进行下一步。

4.2 编译amplfunc.mexw32文件

http://blog.csdn.net/ky25103378/article/details/50496203 博客中详述了如何在VS中编译mex文件的方法。其中几个注意地方:
1. 建立的事DLL工程
2. 打开项目属性配置页，C++ -> 附加包含目录加入MATLAB安装目录下的 \extern\include 路径。
3. 链接器 -> 附加库目录加入MATLAB的 \extern\lib\win32\microsoft 路径。
4. 链接器 -> 输入 -> 附加依赖项输入libmx.lib libeng.lib libmat.lib libmex.lib 这四个lib文件。
5. 输出文件名改为*.mexw32
6. 我的VS是32位的，因此整个过程中都加入的32位source文件
将amplfunc.c加入到工程source文件夹下。点击build，如果一切顺利在解决方案debug目录下出现了amplfunc.mexw32的文件，将此文件拷贝到matlab工作路径下，可进行测试。
（amplfunc.c是针对稠密矩阵，要编译针对稀疏矩阵的需要使用spamfunc.c文件。）

5. 测试

.nl 文件的获得
.nl 文件是AMPL 模型生成的，生成的指令是 write “b+filename”;
如write “bOPFBRN”;会在mod文件所在目录下生成OPFBRN.nl文件，这个文件就可以作为matlab的测试文件了。

在example目录下，利用matlab运行代码

x = 0;
pname = 'OPFBRN'
[x,bl,bu,v,cl,cu] = spamfunc(pname);
[f,c] = spamfunc(x,0);
[g,Jac] = spamfunc(x,1);
W =(spamfunc([1; 1; 1]));
[g,A,B] = evalg(x)
W = evalw(1,[1 1])

如果一切顺利会出现想要的结果。该编译出来的只能在matlab 32位条件下运行，64位会出错。

编译的结果存在了我的CSDN下，如何使用可以参考AMPL Hooking Your Solver to AMPL中的使用说明书。整个编译过程中的文件可以在：http://url.cn/5ywSGKN 下载（好像有效期是40天）。

如有问题欢迎相互讨论。

附：matlab 调用函数

Sample mex function (in MATLAB 5.x format) for getting functions,gradients, and dense Hessians from an AMPL .nl file.  Start with[x,bl,bu,v,cl,cu] = amplfunc('stub')or, for complementarity problems,[x,bl,bu,v,cl,cu,cv] = amplfunc('stub')to read in a problem (discarding the previous problem, if any).The return values are:x = primal initial guessv = dual initial guessbl, bu = lower and upper bounds on xcl, cu = lower and upper bounds on c (the constraint bodies).cv variables complementing constraints:  if cv(i) > 0, thenconstraint i complements x(cv(i)); otherwiseconstraint i is an ordinary constraint.Then[f,c] = amplfunc(x,0)gives the function (f) and constraint bodies (c) at x;[g,Jac] = amplfunc(x,1)gives the gradient g of f, the Jacobian matrix J of c, andW = amplfunc(v)gives the Hessian W of the Lagrangian function L = f + v*c(at the last x at which amplfunc(x,1) was called).After finding optimal values for x and v,amplfunc('solution message',x,v)to write a stub.sol file.Use with MATLAB
It is easy to use AMPL with MATLAB — with the help of a mex file that reads stub.nl files, writes
stub.sol files, and provides function, gradient, and Hessian values. Example file amplfunc.c is source
for an amplfunc.mex that looks at its left- and right-hand sides to determine what it should do and
works as follows:
[x,bl,bu,v,cl,cu] = amplfunc(’stub’)
reads stub.nl and sets
x = primal initial guess,
bl = lower bounds on the primal variables,
bu = upper bounds on the primal variables,
v = dual initial guess (often a vector of zeros),
cl = lower bounds on constraint bodies, and
cu = upper bounds on constraint bodies.
[f,c] = amplfunc(x,0)
sets
f = value of first objective at x and
c = values of constraint bodies at x.
[g,Jac] = amplfunc(x,1)
sets
g = gradient of first objective at x and
Jac = Jacobian matrix of constraints at x.
W = amplfunc(Y)
sets W to the Hessian of the Lagrangian (equation (*) in the section ‘‘Evaluating Nonlinear Functions’’
above) for the first objective at the point x at which the objective and constraint bodies were most recently
evaluated. Finally,
[] = amplfunc(msg,x,v)
calls write_sol(msg,x,v,0) to write the stub.sol file, with
msg = termination message (a string),
x = optimal primal variables, and
v = optimal dual variables.
It is often convenient to use .m files to massage problems to a desired form. To illustrate this, the
examples directory offers the following files (which are simplified forms of files used in joint work with
Michael Overton and Margaret Wright):
• init.m, which expects variable pname to have been assigned a stub (a string value), reads stub.nl,
and puts the problem into the form
minimize f(x)
s. t. c(x) = 0
and d(x) ≥ 0.
For simplicity, the example init.m assumes that the initial x yields d(x) > 0. A more elaborate version
of init.m is required in general.
• evalf.m, which provides [f,c,d] = evalf(x).
October 5, 2000
- 23 -
• evalg.m, which provides [g,A,B] = evalg(x), where A = c’(x) and B = d’(x) are the
Jacobian matrices of c and d.
• evalw.m, which computes the Lagrangian Hessian, W = evalw(y,z), in which y and z are vectors
of Lagrange multipliers for the constraints
c(x) = 0
and
d(x) ≥ 0,
respectively.
• enewt.m, which uses evalf.m, evalg.m and evalw.m in a simple, non-robust nonlinear
interior-point iteration that is meant mainly to illustrate setting up and solving an extended system involving
the constraint Jacobian and Lagrangian Hessian matrices.
• savesol.m, which writes file stub.sol to permit reading a computed solution into an AMPL session.
• hs100.amp, an AMPL model for test problem 100 of Hock and Schittkowski [13].
• hs100.nl, derived from hs100.amp. To solve this problem, start MATLAB and type
pname = ’hs100’;
init
enewt
savesol
Amplfunc.c provides dense Jacobian matrices and Lagrangian Hessians; spamfunc.c is a variant
that provides sparse Jacobian matrices and Lagrangian Hessians. To see an example of using
spamfunc, change all occurrences of ‘‘amplfunc’’ to ‘‘spamfunc’’ in the .m files.