快速使用自定义的非线性求解算法

---

非线性代数系统的通用形式为：

其中， 是变量向量， 是目标函数向量。非线性求解算法的目标就是解算出 。

接下来，我们将通过一个快速示例来演示如何运行您的第一个非线性求解算法。

拷贝文件

若您已经在快速使用自定义的积分算法中运行了自定义积分算法MyEuler，那么mws_user_alg.c文件中实现的非线性求解算法userNonLinearAlg也已经生效，无需再次拷贝该文件。否则，您需要按照快速使用自定义的积分算法的相关步骤拷贝mws_user_alg.c文件。

创建模型

参照创建模型的步骤创建一个新模型，复制以下模型代码：

model Model5
parameter Real p = 1.5;
  Real u;
  Real v;
  Real x;
  Real y;
  Real w = 2 * u;
equation
  u = 0.1 + 2 * sin(3 * time);
  v = 0.5 + cos(5 * time);
  p * x + p * y = 2 * p * u + p * v;
  x ^ 2 - 3 * p * y = 4 * u ^ 2 - 3 * p * v;
end Model5

提示

• 上述示例中的模型名为 Model5，您可以根据需要自定义模型名称；

• 只有当非线性方程系统撕裂后的变量数大于 0 时，仿真时才会调用非线性系统求解算法。如图所示，在上述示例模型翻译后，您可以从输出窗口的建模标签页中看到非线性方程系统撕裂后的变量数为 2。

  

开始仿真

单击仿真标签页上方的仿真按钮，Sysplorer 将调用自定义的非线性求解算法 userNonLinearAlg 来求解模型中的非线性方程系统。至此，您的第一个自定义非线性算法已经成功运行。

提示

与积分算法不同，拷贝完mws_user_alg.c文件之后，其中的自定义非线性系统求解算法将自动生效，无需在界面手动选择。因此，您可能不会明显察觉到这一过程。

自定义非线性求解算法解析

mws_user_alg.c文件内实现了一个经典的牛顿迭代非线性算法，名为 userNonLinearAlg，并将其成功设置到了 Sysplorer 的求解器中。自定义非线性算法一旦设置成功，则系统会自动调用该算法来求解模型中的非线性代数系统。接下来，我们将深入探讨mws_user_alg.c中非线性求解算法实现细节，为您后续扩展自定义非线性求解算法打下坚实的基础。

注意

自定义积分算法是注册到求解器中供用户选择的，而自定义非线性求解算法则是直接设置到求解器中，以替换默认的非线性算法。

常用数据类型

Sysplorer 常用的内部数据类型在各种算法中均一致，详细信息，请参见常用数据类型。完整的 Sysplorer 内部数据类型和宏的定义参见：

%Sysplorer安装目录%\Simulator\Src\mws_common_decl.h

非线性求解算法对象的定义

typedef struct
{
    MwsInteger                  n;
    MwsReal                      tol;
    MwsReal*                     rwork;
    MwsInteger*                 ipiv;
    MwsNLPCallbackFns          cbfns;
    MwsUserContext              userContext;
} UserNLPData;

UserNLPData 各成员的含义如下：

n	求解变量的个数
tol	非线性求解的收敛精度
rwork	存放非线性算法工作数据的数组
ipiv	记录线性方程求解过程中系数矩阵进行的行交换
cbfns	非线性求解算法可直接使用的回调函数集，由求解器提供
userContext	用户数据环境，调用回调函数时需要此参数

创建非线性求解算法的实例

static MwsNLPSolverType UserNLPSolverInstantiate(MwsSize n, const MwsNLPCallbackFns* cbfns, MwsBoolean loggingOn, MwsLocale locale, MwsUserContext userContext)
{
    UserNLPData* userNLPData = cbfns->allocateMemory(userContext, 1, sizeof(UserNLPData));
    if (userNLPData == MWSnullptr)
    {
        return MWSnullptr;
    }

    userNLPData->rwork = cbfns->allocateMemory(userContext, n * (n + 2), sizeof(MwsReal));
    if (userNLPData->rwork == MWSnullptr)
    {
        cbfns->freeMemory(userContext, userNLPData);
        return MWSnullptr;
    }

    userNLPData->ipiv = cbfns->allocateMemory(userContext, n, sizeof(MwsInteger));
    if (userNLPData->ipiv == MWSnullptr)
    {
        cbfns->freeMemory(userContext, userNLPData->rwork);
        cbfns->freeMemory(userContext, userNLPData);
        return MWSnullptr;
    }

    userNLPData->n = n;
    userNLPData->cbfns = *cbfns;
    userNLPData->userContext = userContext;
    return userNLPData;
}

与创建 MyEuler 算法的实例类似，详细信息，请参见创建积分算法的实例，创建 userNonLinearAlg 实例的过程也包括为成员分配内存以及回调函数。

获取非线性算法的特性

static MwsInteger UserNLPSolverGetFeature(MwsNLPSolverType inst, MwsNLPSolverFeature* solverFeature)
{
    MwsInteger rst = 0;
    solverFeature->columnMajorMatrix = mwsTrue;
    solverFeature->gradientNeeded = mwsTrue;
    return rst;
}

UserNLPSolverGetFeature 函数用于获取非线性求解算法的特性。与前述的积分算法 MyEuler 相比，非线性求解算法的特性要简单得多，仅包含两个特性。通常情况下，columnMajorMatrix 设置为 mwsTrue。gradientNeeded 根据算法实际情况设置。由于 userNonLinearAlg 使用的是牛顿迭代法，所以 gradientNeeded 设置为 mwsTrue。

columnMajorMatrix	矩阵是否为列主排列
gradientNeeded	是否需要梯度信息

设置非线性算法的参数

static MwsInteger UserNLPSolverSetExperiment(MwsNLPSolverType inst, MwsReal tolerance, MwsBoolean maximumTimesDefined, MwsInteger maximumTimes, void* reserve)
{
    MwsInteger rst = 0;
    UserNLPData* userNLPData = (UserNLPData*)(inst);
    userNLPData->tol = tolerance;
    return rst;
}

UserNLPSolverSetExperiment 用来给非线性算法实例 inst 设置求解配置，可以控制非线性求解的过程。例如，参数 tolerance 可以指定求解的收敛精度，参数 maximumTimes 指定最大迭代次数。

有关 UserNLPSolverSetExperiment 的详细接口说明，请参见设置非线性算法的参数。

求解非线性代数系统

static MwsInteger UserNLPSolve(MwsNLPSolverType inst, const MwsChar* const xname[], const MwsReal xnominal[], MwsReal x[], void* reserve)
{
    MwsInteger rst = 0;
    UserNLPData* userNLPData = (UserNLPData*)(inst);
    MwsInteger n = userNLPData->n;
    MwsInteger nrhs = 1;
    MwsInteger info = 0;
    MwsInteger i = 0;
    MwsInteger maxIter = 50;
    MwsInteger iterNum = 0;
    MwsReal tol = 1e-5;
    MwsReal fvalFac = 0.0;
    MwsReal* fval = userNLPData->rwork;
    MwsReal* dval = fval + n;
    MwsReal* jacVal = dval + n;


    rst = userNLPData->cbfns.fncEvaluation(userNLPData->userContext, x, fval);
    if (rst != 0)
    {
        return rst;
    }

    do
    {
        rst = userNLPData->cbfns.jacEvaluation(userNLPData->userContext, x, jacVal);
        if (rst != 0)
        {
            break;
        }

        memcpy(dval, fval, n * sizeof(MwsReal));
        for (i = 0; i < n; i++)
        {
            dval[i] *= -1;
        }

        rst = dgesv_(&n, &nrhs, jacVal, &n, userNLPData->ipiv, dval, &n, &info);
        if (rst != 0)
        {
            break;
        }

        for (i = 0; i < n; i++)
        {
            x[i] += dval[i];
        }

        rst = userNLPData->cbfns.fncEvaluation(userNLPData->userContext, x, fval);
        if (rst != 0)
        {
            break;
        }

        fvalFac = 0.0;
        for (i = 0; i < n; i++)
        {
            fvalFac += fval[i] * fval[i];
        }

        iterNum++;

    } while (sqrt(fvalFac) > tol || iterNum <= maxIter);

    return rst;
}

UserNLPSolve 各参数的意义如下：

inst	非线性算法实例
xname	求解变量的名字向量（主要用于日志输出）
xnominal	求解变量的标称值（由modelica规范定义）
x	求解变量值（输入时为初始值，输出时为结果值）
reverse	保留参数

UserNLPSolve 实现了牛顿迭代法的核心逻辑，具体解释如下：

1. 首先调用 fncEvaluation 计算 的函数值然后开始 do 循环；
2. 在循环体内，先调用 jacEvaluation 计算 的雅可比矩阵 ；
3. 接着调用外部函数 dgesv_求解线性方程 ，获取求解变量的增量 ；dgesv_函数的说明，请参见求解线性代数系统；
4. 然后利用该增量更新求解变量 ；
5. 最后检查 是否满足精度要求，如果满足则终止迭代，否则进入下一轮循环继续迭代。

释放算法实例

static void UserNLPSolverFreeInstance(MwsNLPSolverType inst)
{
    UserNLPData* userNLPData = (UserNLPData*)(inst);
    if (userNLPData != MWSnullptr)
    {
        userNLPData->cbfns.freeMemory(userNLPData->userContext, userNLPData->rwork);
        userNLPData->cbfns.freeMemory(userNLPData->userContext, userNLPData->ipiv);
        userNLPData->cbfns.freeMemory(userNLPData->userContext, userNLPData);
        userNLPData = MWSnullptr;
    }
}

UserLPSolverFreeInstance 函数用于释放算法实例的内存。

设置非线性算法

MwsStatus MwsSetExternalNonlinearSolver(void* inst)
{
    /* Set user defined Nonlinear solver by calling following function. */
    /* MwsStatus MwsSetUserNonlinearSolver(void* inst, const MwsString solverName, const MwsNLPSolverFcns* nlpSolverFcns) */
    /* you must provide valid callback functions in nlpSolverFcns */
    /* you can only register one nonlinear solver, the last one will be effective if more than one registered */

    MwsStatus rst = mwsOK;
    MwsNLPSolverFcns nlpSolverFcns;
    MwsString solverName = "userNonLinearAlg";
    nlpSolverFcns.solverInstantiate = UserNLPSolverInstantiate;
    nlpSolverFcns.solverGetFeature = UserNLPSolverGetFeature;
    nlpSolverFcns.solverSetExperiment = UserNLPSolverSetExperiment;
    nlpSolverFcns.solverSetDebugLogging = UserNLPSolverSetDebugLogging;
    nlpSolverFcns.solverSolve = UserNLPSolve;
    nlpSolverFcns.solverFreeInstance = UserNLPSolverFreeInstance;
    rst = MwsSetUserNonlinearSolver(inst, solverName, &nlpSolverFcns);
    return rst;
}

MwsSetExternalNonLinearSolver 函数用于将 Sysplorer 求解器中的非线性求解算法设置为自定义的非线性求解算法。只有在算法成功设置后，Sysplorer 才能正确识别并执行用户自定义的算法。设置非线性求解算法的主要步骤与注册积分算法的步骤极为相似，因此不再赘述。 有关MwsSetExternalNonLinearSolver函数的使用说明，请参见用户手册的设置函数的模板。