makeFuncFile函数
|
5
|
0
|
@DynamicModel

885 字
|
4 分钟

这段代码定义了一个函数 makeFuncFile，用于将 DynamicModel 对象转换为一个 MATLAB 函数文件，以便快速求解模型。以下是代码的详细解释：

1. 函数功能

makeFuncFile 函数的主要功能是：

  • 将 DynamicModel 对象转换为一个 MATLAB 函数文件。

  • 该函数文件可以直接用于求解模型，而不需要使用 DynamicModel 对象和 DynamicElement 等结构。

  • 生成的函数文件包含模型的微分方程、事件处理逻辑以及输入、控制、参数等的处理。

2. 输入参数

  • obj：DynamicModel 对象，表示动态模型。

  • filename：生成的 MATLAB 函数文件的路径和名称（包括 .m 后缀）。

  • funcName：生成的 MATLAB 函数的名称（通常与文件名相同，不包括路径）。

3. 函数逻辑

3.1 获取字段名称

matlab
复制
[stateNames, auxNames, ctrlNames, paramNames, inputNames] = getFieldNames(obj);

  • 使用 getFieldNames 函数获取 DynamicModel 对象中所有状态变量、辅助变量、控制变量、参数和输入变量的名称。

3.2 创建函数文件

matlab
复制
fid = fopen(filename, 'w');

  • 打开文件 filename 以写入内容。

3.3 写入函数头

matlab
复制
fprintf(fid, ['function ' funcName '(obj, solver, options)\n\n']);

  • 写入函数头，定义函数名和输入参数。

3.4 处理输入参数

matlab
复制
fprintf(fid, ['if ~exist(''options'',''var'')\n' ...
    '\t options = [];\n' ...
    'end\n\n']);

  • 如果未提供 options 参数，则将其设置为空。

3.5 获取字段名称

matlab
复制
fprintf(fid, ['[stateNames, auxNames, ctrlNames, ...\n'...
    '\t paramNames, inputNames] = getFieldNames(obj);\n']);

  • 在生成的函数中调用 getFieldNames 函数，获取字段名称。

3.6 设置时间变量

matlab
复制
fprintf(fid, ['\n%% Set time variable\n' ...
    'time = obj.t.val;\n']);

  • 设置时间变量 time。

3.7 转换输入变量为矩阵

matlab
复制
fprintf(fid, '\n%% Convert inputs to matrix\n');
fprintf(fid, 'd = obj.d.(inputNames{1}).val(:,1);\n');
fprintf(fid, 'for k=1:length(inputNames)\n');
fprintf(fid, '\t d = [d obj.d.(inputNames{k}).val(:,2)];\n');
fprintf(fid, 'end\n');

  • 将输入变量转换为矩阵形式。

3.8 转换预定义控制变量为矩阵

matlab
复制
fprintf(fid, '\n%% Convert predefined controls to matrix\n');
fprintf(fid, 'u = [];\n');
fprintf(fid, 'uTime = [];\n');
fprintf(fid, 'for k=1:length(ctrlNames)\n');
fprintf(fid, ['\t if isnumeric(obj.u.(ctrlNames{k}).val) && ~isscalar(obj.u.(ctrlNames{k}).val) && ~isempty(obj.u.(ctrlNames{k}).val)\n' ...
    '\t\t uTime = obj.u.(ctrlNames{k}).val(:,1);\n' ...
    '\t\t if sum(isnan(u(:))) == size(u,1)*size(u,2) \n' ...
           '\t\t\t u = nan(length(obj.u.(ctrlNames{k}).val(:,2)), length(u));\n' ...
     '\t\t end\n' ...
    '\t\t u = [u obj.u.(ctrlNames{k}).val(:,2)];\n' ...
    '\t else\n' ...
        '\t\t if isempty(u)\n'...
            '\t\t\t u = nan(1,1);\n' ...
        '\t\t else\n' ...
            '\t\t\t u = [u nan(length(u(:,1)),1)];\n' ...
        '\t\t end\n' ...
    '\t end\n']);
fprintf(fid, 'end\n');
fprintf(fid, ['if ~isempty(uTime)\n' ...
    '\t u = [uTime u];\n' ...
    'end\n']);

  • 将预定义的控制变量转换为矩阵形式。

3.9 创建参数数组

matlab
复制
fprintf(fid, '\n%% Create array of parameters\n');
for k=1:length(paramNames)
    if defIsLabel(obj.p.(paramNames{k})) % Parameter has no definition
        val = obj.p.(paramNames{k}).val;
    else % Parameter is based on other parameters
        val = obj.p.(paramNames{k}).def([],[],[],[],p);
        obj.p.(paramNames{k}).val = val; 
    end
    p.(paramNames{k}) = val;
    fprintf(fid, ['p(' num2str(k) ') = ' ...
        num2str(val) ';\n']);
end

  • 创建参数数组 p。

3.10 运行仿真

matlab
复制
fprintf(fid, '\n%% Run simulation\n');
if isempty(obj.e) % no events defined
    fprintf(fid, '[t,x] = feval(solver, @(t,x) ode(x, t, d, u, p), time, getInitialStates(obj), options);\n');
    fprintf(fid, 'setSolution(obj, t, x);\n\n');
    fprintf(fid, 'end\n');
else % include an event listener in the ODE solver
    fprintf(fid, 'options = odeset(options, ''Events'', @events);\n');
    fprintf(fid, 'tOut = [];\n');
    fprintf(fid, 'xOut = [];\n');
    fprintf(fid, 'tStart = time(1);\n');
    fprintf(fid, 'tFinal = time(2);\n');
    fprintf(fid, 'xInit =  getInitialStates(obj);\n');
    fprintf(fid, 'refine = odeget(options,''Refine'');\n');
    fprintf(fid, 'if isempty(refine)\n'); 
        fprintf(fid, '\trefine = 1;\n');
    fprintf(fid, 'end\n\n');
    
    % Default for MaxStep is 0.1*abs(tFinal-tStart)
    fprintf(fid, '%% Default for MaxStep is 0.1*abs(tFinal-tStart)\n');
    fprintf(fid, 'if isempty(odeget(options,''MaxStep'')) %% no MaxStep has been defined\n');
        fprintf(fid, '\t options = odeset(options, ''MaxStep'', 0.1*abs(tFinal-tStart));\n');
    fprintf(fid, 'end\n\n'); % function
    
    for k=1:length(obj.e)
        varString = '';
        valString = '';
        for n=1:length(obj.e(k).resetVars)
            label = obj.e(k).resetVars(n).label;
            label = label(3:end);
            varString = [varString ' ' num2str(find(strcmp(stateNames,label)))];
            valString = [valString ' ' num2str(obj.e(k).resetVals(n))];
        end
        fprintf(fid, ['eventVars{%d} = [' varString '];\n'], k);
        fprintf(fid, ['eventVals{%d} = [' valString '];\n'], k);
    end
            
    % Iteratively solve until event is reached. Then reset values
    % according to event definitions
    fprintf(fid, '\n\n%% Iteratively solve until event is reached. Then reset values\n');
    fprintf(fid, '%% according to event definitions\n');
    
    fprintf(fid, 'while tStart < tFinal\n');
        fprintf(fid, '\t [t,x,~,~,eventNum] = feval(solver, @(t,x) ode(x, t, d, u, p), [tStart tFinal], xInit, options);\n\n');

        % Accumulate output
        fprintf(fid, '\t %% Accumulate output\n');
        fprintf(fid, '\t nt = length(t);\n');
        fprintf(fid, '\t tOut = [tOut; t(2:nt)];\n');
        fprintf(fid, '\t xOut = [xOut; x(2:nt,:)];\n\n');
        fprintf(fid, '\t tStart = t(nt);\n');

        fprintf(fid, '\t if ~isempty(eventNum) %% simulation not finished\n');
            % Set the new initial conditions
            fprintf(fid, '\t\t %% Set the new initial conditions\n');
            fprintf(fid, '\t\t xInit = x(nt,:);\n');
            fprintf(fid, '\t\t xInit(eventVars{eventNum(1)}) = eventVals{eventNum(1)};\n\n');

            % Guess of a first timestep is length of the last valid timestep
            fprintf(fid, '\t\t %% Guess of a first timestep is length of the last valid timestep\n');
            fprintf(fid, '\t\t options = odeset(options, ''InitialStep'', t(nt)-t(nt-refine));\n');
        fprintf(fid, '\t end\n'); % if ~isempty(eventNum)
    fprintf(fid, 'end\n\n'); % while
    fprintf(fid, 'setSolution(obj, tOut, xOut);\n\n');

    fprintf(fid, 'end\n\n'); % function

    % Events for ODE solver
    fprintf(fid, '%% Events for ODE solver\n');
    fprintf(fid, 'function [value,isterminal,direction] = events(~,x)\n');
    
        valString = '';
        dirString = '';
            for k=1:length(obj.e)
                label = obj.e(k).condition.label;
                label = label(3:end);
                valString = [valString ' x(' num2str(find(strcmp(stateNames,label))) ')'];
                dirString = [dirString ' ' num2str(obj.e(k).direction)];
            end
        fprintf(fid, ['\t value = [' valString '];\n']);
        fprintf(fid, ['\t direction = [' dirString '];\n']);
        fprintf(fid, '\t isterminal = ones(1,%d);\n', length(obj.e));
    fprintf(fid, 'end\n'); % events
    
end

  • 如果模型中没有定义事件，则直接运行仿真。

  • 如果模型中定义了事件，则使用事件处理逻辑迭代求解。

3.11 定义 ODE 函数

matlab
复制
fprintf(fid, '\n%% Function for the ODE solver\n');
fprintf(fid, 'function dx = ode(x, t, d, u, p)\n\n');

  • 定义 ODE 函数 ode，用于计算状态变量的导数。

3.12 处理输入和控制变量

matlab
复制
fprintf(fid, ['\t %% sample inputs at time t\n' ...
    '\t if t < d(1,1)\n' ...
        '\t\t dSample = d(1,2:end);\n' ...
    '\t elseif t > d(end,1)\n' ...
        '\t\t dSample = d(end,2:end);\n' ...
    '\t else \n' ...
    '\t\t dSample = nan(size(d,2)-1,1);\n' ...
    '\t\t for k=1:(size(d,2)-1)\n' ...
        '\t\t\t dSample(k) = interp1(d(:,1),d(:,k+1),t);\n'...
    '\t\t end\n' ...
    '\t end\n' ...
    '\t d = dSample;\n']);

  • 在时间 t 处对输入变量进行插值。

matlab
复制
fprintf(fid, ['\n\t %% sample predefined controls at time t\n' ...
    '\t if t < u(1,1)\n' ...
        '\t\t uSample = u(1,2:end);\n' ...
    '\t elseif t > u(end,1)\n' ...
        '\t\t uSample = u(end,2:end);\n' ...
    '\t else \n' ...
    '\t\t uSample = nan(size(u,2)-1,1);\n' ...
    '\t\t for k=1:(size(u,2)-1)\n' ...
        '\t\t\t if ~isnan(u(1,k+1))\n' ...
            '\t\t\t\t uSample(k) = interp1(u(:,1),u(:,k+1),t);\n'...
        '\t\t\t end\n' ...
    '\t\t end\n' ...
    '\t end\n' ...
    '\t u = uSample;\n']);

  • 在时间 t 处对预定义的控制变量进行插值。

3.13 计算辅助变量

matlab
复制
fprintf(fid, '\n\t%% Create struct of auxiliary states\n');
for k=1:length(auxNames)
    auxDef = remStructs(obj, obj.a.(auxNames{k}).def);       
    fprintf(fid, ['\t a.' auxNames{k} ' = ' ...
        auxDef ';\n']);
end

  • 计算辅助变量的值。

3.14 计算状态变量的导数

matlab
复制
fprintf(fid, '\n\t%% Create array of derivatives \n');
for k=1:length(stateNames)
    odeDef = remStructs(obj, obj.x.(stateNames{k}).def);       
    fprintf(fid, ['\t dx(' num2str(k) ') = ' ...
        odeDef ';\n']);
end

  • 计算状态变量的导数。

3.15 结束 ODE 函数

matlab
复制
fprintf(fid, '\n\t dx = dx'';\n');
fprintf(fid, 'end\n');

  • 结束 ODE 函数。

3.16 关闭文件

matlab
复制
fclose(fid);

  • 关闭文件。

4. 辅助函数 remStructs

remStructs 函数用于替换字符串中的参数、输入和控制变量名称，将其替换为数组值。

matlab
复制
function newStr = remStructs(obj, str)
    if isa(str, 'function_handle')
        str = func2str(str);
        str = str(13:end);
    end

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

    newStr = str;
    % replace all parameter names with their array value
    for n=1:length(paramNames)
        paramLengths(n) = length(paramNames{n});
    end
    [~, paramByLength] = sort(paramLengths,'descend');
    for n=paramByLength
        newStr = strrep(newStr, ['p.' paramNames{n}], ['p(' num2str(n) ')']);
    end

    % replace all input names with their array value
    for n=1:length(inputNames)
        inputLengths(n) = length(inputNames{n});
    end
    [~, inputByLength] = sort(inputLengths,'descend');
    for n=inputByLength
        newStr = strrep(newStr, ['d.' inputNames{n}], ['d(' num2str(n) ')']);
    end

    % replace all predefined control names with their array value
    if ~isempty(ctrlNames)
        for n=1:length(ctrlNames)
            ctrlLengths(n) = length(ctrlNames{n});
        end
        [~, ctrlByLength] = sort(ctrlLengths,'descend');
        for n=ctrlByLength
            newStr = strrep(newStr, ['u.' ctrlNames{n}], ['u(' num2str(n) ')']);
        end
    end
    
    % replace all states names with their array value
    for n=1:length(stateNames)
        stateLengths(n) = length(stateNames{n});
    end
    [~, stateByLength] = sort(stateLengths,'descend');
    for n=stateByLength
        newStr = strrep(newStr, ['x.' stateNames{n}], ['x(' num2str(n) ')']);
    end
end

5. 总结

makeFuncFile 函数的作用是将 DynamicModel 对象转换为一个 MATLAB 函数文件，以便快速求解模型。生成的函数文件包含模型的微分方程、事件处理逻辑以及输入、控制、参数等的处理。该函数适用于需要将动态模型转换为独立 MATLAB 函数的场景，例如提高求解效率或与其他工具集成。

暂无评论

发送评论 编辑评论 


				

			

						

				

					

						

							

								
							

							
						

					

				

				

					

						

							

								
							

							
						

					

				

				

					

						

							

								
							

							
							
													

					

				

			

			
			

				

					
				

			

				

											

							
							
						

																				
					
											
						

	

		
|´・ω・)ノ
ヾ(≧∇≦*)ゝ
(☆ω☆)
（╯‵□′）╯︵┴─┴
￣﹃￣
(/ω＼)
∠( ᐛ 」∠)＿
(๑•̀ㅁ•́ฅ)
→_→
୧(๑•̀⌄•́๑)૭
٩(ˊᗜˋ*)و
(ノ°ο°)ノ
(´இ皿இ｀)
⌇●﹏●⌇
(ฅ´ω`ฅ)
(╯°A°)╯︵○○○
φ(￣∇￣o)
ヾ(´･ ･｀｡)ノ"
( ง ᵒ̌皿ᵒ̌)ง⁼³₌₃
(ó﹏ò｡)
Σ(っ °Д °;)っ
( ,,´･ω･)ﾉ"(´っω･｀｡)
╮(╯▽╰)╭
o(*////▽////*)q
＞﹏＜
( ๑´•ω•) "(ㆆᴗㆆ)
😂
😀
😅
😊
🙂
🙃
😌
😍
😘
😜
😝
😏
😒
🙄
😳
😡
😔
😫
😱
😭
💩
👻
🙌
🖕
👍
👫
👬
👭
🌚
🌝
🙈
💊
😶
🙏
🍦
🍉
😣
Source: github.com/k4yt3x/flowerhd
	

	

		
颜文字
Emoji
小恐龙
花!
	


									

			

			
			
		

	






上一篇
下一篇