odeSolve函数 – WordPress

这段代码定义了一个名为 odeSolve 的 MATLAB 函数，用于对状态空间模型（obj）进行数值求解。与 ode15s 类似，它通过调用 MATLAB 的 ODE 求解器（如 ode15s 或其他求解器）来求解状态空间模型的微分方程。以下是代码的详细解释：

函数签名

function [t, x] = odeSolve(obj, func, options)

输入参数:
- obj: 状态空间模型对象，包含状态、辅助状态、控制、参数和输入的定义。
- func: 指定的 ODE 求解器函数（如 ode15s、ode45 等）。
- options: 可选的 ODE 求解器选项（如相对误差、绝对误差等）。
输出参数:
- t: 时间向量，表示仿真时间点。
- x: 状态变量的值矩阵，每一列对应一个状态变量在不同时间点的值。

代码逻辑

1. 检查输入参数

if ~exist('options','var')
    options = [];
end

如果未提供 options 参数，则将其初始化为空数组。

2. 获取字段名称

[stateNames, auxNames, ctrlNames, ...
    paramNames, inputNames] = getFieldNames(obj);

从状态空间模型对象 obj 中获取状态、辅助状态、控制、参数和输入的名称列表。

3. 提取辅助状态函数

for n=1:length(auxNames)
    a.(auxNames{n}) = obj.a.(auxNames{n}).func;
end

从 obj 中提取辅助状态的函数句柄，并存储在结构体 a 中。

4. 提取控制函数或值

for n=1:length(ctrlNames)
    if strcmp(obj.u.(ctrlNames{n}).def, ['u.' ctrlNames{n}]) 
        % the def for the control is the control name, control acts as input
        u.(ctrlNames{n}) = obj.u.(ctrlNames{n}).val;
    else % control is rule based, acts as an auxiliary state
        defExpand(obj, obj.u.(ctrlNames{n}));
        u.(ctrlNames{n}) = obj.u.(ctrlNames{n}).func;
    end
end

遍历控制名称列表：
- 如果控制的定义是控制名称本身（即控制作为输入），则直接提取其值。
- 如果控制是基于规则的（即控制作为辅助状态），则通过 defExpand 函数扩展定义，并提取其函数句柄。

5. 提取输入数据

for n=1:length(inputNames)
    d.(inputNames{n}) = obj.d.(inputNames{n}).val;
end

从 obj 中提取输入数据，并存储在结构体 d 中。

6. 提取参数值

for n=1:length(paramNames)
    p.(paramNames{n}) = obj.p.(paramNames{n}).val;
end

从 obj 中提取参数值，并存储在结构体 p 中。

7. 提取状态变量的导数函数

for n=1:length(stateNames)
    xOde.(stateNames{n}) = obj.x.(stateNames{n}).func;
end

从 obj 中提取状态变量的导数函数句柄，并存储在结构体 xOde 中。

8. 调用 ODE 求解器

[t,x] = feval(func, @(t,x) getOdes(t, x, a, u, d, p, xOde, stateNames), obj.t.val, getInitialStates(obj), options);

使用 feval 调用指定的 ODE 求解器 func，求解状态空间模型的微分方程。
- @(t,x) getOdes(t, x, a, u, d, p, xOde, stateNames): 匿名函数，用于计算状态变量的导数。
- obj.t.val: 时间范围（仿真时间点）。
- getInitialStates(obj): 获取状态变量的初始值。
- options: ODE 求解器的选项。

9. 设置仿真结果到状态变量

for n=1:length(stateNames)
    obj.x.(stateNames{n}).val = [t x(:,n)];
    xStruct.(stateNames{n}) = x(:,n);
end
x = xStruct;

将仿真结果（时间 t 和状态变量值 x）保存到 obj 的状态变量中。
同时将状态变量值存储在结构体 xStruct 中，并将其赋值给输出变量 x。

10. 插值输入数据

for n=1:length(inputNames)
    d.(inputNames{n}) = interp1(obj.d.(inputNames{n}).val(:,1),...
        obj.d.(inputNames{n}).val(:,2),t);    
end

根据仿真时间 t 对输入数据进行插值，得到每个时间点的输入值，并存储在结构体 d 中。

11. 计算并设置控制值

for n=1:length(ctrlNames)
    if isscalar(obj.u.(ctrlNames{n}).val) || isempty(obj.u.(ctrlNames{n}).val)
        % control was not predefined, need to calculate
        try
            defExpand(obj, obj.u.(ctrlNames{n}));
            if ~exist('a','var')
                a = [];
            end
            if ~exist('u','var')
                u = [];
            end
            u.(ctrlNames{n}) = obj.u.(ctrlNames{n}).func(x, a, u, d, p);
            if isscalar(u.(ctrlNames{n})) % the definition does not depend on t, it's constant
                u.(ctrlNames{n}) = u.(ctrlNames{n})*ones(size(t));
            end
            obj.u.(ctrlNames{n}).val = [t u.(ctrlNames{n})];
        catch err
            msg = sprintf('%s \n\nFailed to evaluate the definition of DynamicElement u.%s (k=%d): \n\t''%s''',...
                err.message, ctrlNames{n}, n, obj.u.(ctrlNames{n}).def);
            id = 'MATLAB:StateSpaceModel:eval';
            error(id,msg);
        end
    else % interpolate u in timepoints t
        u.(ctrlNames{n}) = interp1(obj.u.(ctrlNames{n}).val(:,1),...
            obj.u.(ctrlNames{n}).val(:,2),t); 
    end
end

遍历控制名称列表，计算或插值控制值：
- 如果控制值是标量或未定义，则通过 defExpand 函数扩展定义，并调用其函数句柄计算控制值。
- 如果控制值是时间相关的，则根据仿真时间 t 进行插值。
将计算结果保存到 obj 的控制变量中。

12. 计算并设置辅助状态值

for n=1:length(auxNames)
    try
        if ~exist('a','var')
            a = [];
        end
        if ~exist('u','var')
            u = [];
        end
        a.(auxNames{n}) = obj.a.(auxNames{n}).func(x, a, u, d, p);
        if isscalar(a.(auxNames{n})) % the definition does not depend on t, it's constant
            a.(auxNames{n}) = a.(auxNames{n})*ones(size(t));
        end
        obj.a.(auxNames{n}).val = [t a.(auxNames{n})];
    catch err
        msg = sprintf('%s \n\nFailed to evaluate the definition of DynamicElement a.%s (n=%d): \n\t''%s''',...
            err.message, auxNames{n}, n, obj.a.(auxNames{n}).def);
        id = 'MATLAB:StateSpaceModel:eval';
        error(id,msg);
    end
end

遍历辅助状态名称列表，通过函数句柄计算辅助状态的值。
如果辅助状态的定义与时间无关，则将其扩展为与时间向量 t 相同长度的常量向量。
将计算结果保存到 obj 的辅助状态变量中。

总结

odeSolve 函数的主要作用是使用指定的 ODE 求解器对状态空间模型进行数值求解，并将仿真结果保存到模型对象中。具体功能包括：

提取模型中的辅助状态、控制、输入、参数和状态变量的导数函数。
调用指定的 ODE 求解器进行仿真。
将仿真结果保存到状态变量中。
计算或插值控制值和辅助状态值，并将其保存到模型对象中。

该函数在状态空间模型的仿真和分析中非常有用，特别是在需要动态计算控制值和辅助状态值时。