Use of User-defined Language Specifications I – Nonlinear equation system with Matlab

Goal

In this example, we will set up a UDLS for Matlab. This UDLS shall generate code to solve a nonlinear equation system containing function calls. The UDLS is available in MOSAICmodeling with ID 183144.

Workflow for setting up a UDLS

A UDLS is set up in the “Miscellaneous” section of MOSAICmodeling in the tab Language Specification. Should this tab be invisible to you and you would like to use it, please request this feature for your account.

Global language settings

We begin in the tab Global Language Settings by adding an informative description. As programming languages can change their syntax over time, you should also mention for which version this UDLS is tested. Then, perform the following steps:

1. Activate Direct Functions Supported and Case Sensitive in the Language Settings. The former will allow us to include direct function calls, the latter ensures a difference between capital and small letters
2. Activate Force Floating Point in the Number settings
3. Select “NLE” as Equation System Type
4. Select “Simulation” in the box Opt or Nopt

The final specifications are shown in Figure 1. The button in the bottom right corner next to Minimum Specifications Complete should be green now.

Operations

Move to the Operations tab. We need to specify in which way basic operations are expressed in the programming language, i.e., Matlab. Therefore,

1. Click on Add and add a “+” as Infix, click on the green arrow to see the output
2. Click on Subtract and add a “-” as Infix, click on the green arrow to see the output
3. Click on Multiply and add “(” as Prefix, “)*(” as Infix, and “)” as Suffix, click on the green arrow to see the output. The additional brackets ensure that any addition or subtraction terms in LEFT or RIGHT are correctly enclosed
4. Click on Divide and add “(” as Prefix, “)/(” as Infix, and “)” as Suffix, click on the green arrow to see the output
5. Click on Power and add “(” as Prefix, “)^(” as Infix, and “)” as Suffix, click on the green arrow to see the output
6. Click on R-root and add “(” as Prefix, “)^(1.0/(” as Infix, and “))” as Suffix, click on the green arrow to see the output. Alternatively, you could implement the Matlab function nthroot(X,N). Attention: Make sure that you selected Swap LR because MOSAICmodeling uses the function \sqrt[N]{X} internally, i.e., the first argument (LEFT) is the root’s exponent. This is not shown in Figure 2
7. Click on Exp and add “exp(” as Prefix and “)” as Suffix, click on the green arrow to see the output
8. Click on Log and add “log(” as Prefix and “)” as Suffix, click on the green arrow to see the output
9. Click on Sinus and add “sin(” as Prefix and “)” as Suffix, click on the green arrow to see the output
10. Click on Cosinus and add “cos(” as Prefix and “)” as Suffix, click on the green arrow to see the output
11. Click on Parentheses and add “(” as Prefix and “)” as Suffix, click on the green arrow to see the output

The remaining operations are not of interest for simulation code and can therefore be neglected. The final specifications are shown in Figure 2.

Variable namings

Go to the tab Variable Namings and do the following:

1. Select Use Variable Naming as Strategy
2. Click on Namespace to base name and add a “_” as Infix, click on the green arrow to see the output. You see that there is now a “_” between std and P
3. Add the same “_” for Leading superscript, Leading subscript, and Leading indexname. After every change, click on the green arrow to see what has changed
4. In this example, we leave the Infix for Index name to index value empty

The final specifications are shown in Figure 3.

Code generation properties

Next, we skip the tab External Functions as we do not supply additional interfaces, and go directly to the tab Code Generation Properties. Now, we want to add a couple of basic solver options for Matlab’s fsolve:

1. Add the the Property MaxIter with Handle _MAXITER_, Type INTEGER, and Default 1000
2. Add the the Property TolFun with Handle _TOLFUN_, Type FLOAT, and Default 1e-10
3. Add the the Property DisplayIter with Handle _DISPLAY_, Type STRING, and Default ‘iter’
4. Add the the Property MaxFunEvals with Handle _MAXFUNEVALS_, Type INTEGER, and Default 100000
5. Add the the Property File Name with Handle _FILENAME_, Type STRING, and Default MATLAB_NLE_export
6. Click on the row containing the name DisplayIter and click on “+ Option” on the right-hand side. Enter ‘iter’ and confirm. Repeat this for ‘notify’, ‘off’, ‘final’, ‘iter-detailed’, and ‘final-detailed’. This allows the user to select a preferred choice within MOSAICmodeling when a simulation shall be exported

The final specifications are shown in Figure 4.

Code elements

This tab is fairly complicated as it contains the structure for the generated code (see Figure 5). We will try to be as concise as possible in our description. However, you can of course also directly open the UDLS with the ID given above.

We begin in the in the block _FULL_PROBLEM_SOLVING_CODE_. This block shall basically contain the overall structure of the generated code. Here, we now add the following lines:

_Help_Header_

_Initialize_And_Solve_

_Model_Equations_

_DIR_FUN_IMPL_LOOP_

_Display_Results_

These blocks will execute the following tasks:

• _Help_Header_: generation of documentation code at the beginning
• _Initialize_And_Solve_: variable initialization and call of fsolve
• _Model_Equations_: model equations
• _DIR_FUN_IMPL_LOOP_: loop writing the code for the functions that are part of the model
• _Display_Results_: will show the results for the calculated variables

We will now go through these blocks in sequential order. Note that when you really develop a UDLS, this process might be much more iterative than in this example here.

Help header

The block _Help_Header_ does not exist so far.  Hence, we create it by highlighting the _FULL_PROBLEM_SOLVING_CODE_ block and clicking on “+ Node” at the bottom left corner. Define the handle _Help_Header_ in the popup window and click on Submit”. You can now move the created block upwards or downwards using the two arrows.

Highlight the new block and enter the following code:

% Matlab code to solve the nonlinear equation system below
% Code generated with MOSAICmodeling

Hence, this block will only add a couple of comments at the beginning. You can make this much more complicated and add the whole notation of your system, but we will keep it simple here.

Initialize and solve

The next block is responsible for variable and parameter specification as well as the call of fsolve. First, we need to create the block again in the same way as above. Then, we highlight it and add the following code:

function InitializeAndSolveEquationSystem()

	% load variable start values
_Set_Start_VARS_LOOP_
	% load parameters
_Set_PARAMS_LOOP_
	options = optimoptions('fsolve','MaxIter',_MAXITER_,'MaxFunEvals',_MAXFUNEVALS_,'TolFun',_TOLFUN_,'Display',_DISPLAY_);
	[SOL,RES,EXITFLAG] = fsolve (@( X_ITER ) modelEquations(X_ITER,PARAMS),X_INIT,options);

	displayResults(SOL);
	fprintf('\n')
	fprintf('Residuals: \n')
	disp(RES')
	fprintf('Terminated with exitflag %i \n',EXITFLAG)

end

This code creates a function in Matlab that begins with a looping setting the initial guesses for the variables. Then, the parameter values are set within a loop. Afterwards, the solver options are set based on the handles that were created in the tab Code Generation Properties. Then, we have the call to fsolve. Finally, we display the results, the residuals, and the exitflag of the solver. As you can see, this is basically Matlab code in combination with a couple of handles and blocks of the UDLS.

As you can see, the block requires the loops _Set_Start_VARS_LOOP_ and _Set_PARAMS_LOOP_. These loops are now added as subnodes of the node _Initialize_And_Solve_. Add the first node as “Variables” block and the second node as “Parameters” block. After you have added both nodes, highlight the first one.

One the right-hand side, you can see all available subhandles, e.g., _CTRVALUE_, which is the index of the variables. Add the code below to this block:

X_INIT(_CTRVALUE_) = _VAR_VALUE_;   % _VAR_NAME_

This code will be executed as loop and will write the variable X_INIT with an increasing index. The respective enry of X_INIT will receive the value _VAR_VALUE_.

Then, go to the PARAMS loop and the following line:

PARAMS(_CTRVALUE_) = _PAR_VALUE_;	% _PAR_NAME_

This will do the same for the parameter values. At every point of your code development, you can go to the tab Overall Results and load any simulation you have access to. Then, you can go back to the tab Code Elements and click on “Code” at the bottom right corner of the editor. This generates the code of the highlighted block according to your current specifications.

Model equations

We now move on to formulating the model equations. Therefore, create the block _Model_Equations_ and add the following code:

function[Y] = modelEquations(X,PARAMS)

	% read out variables  
_Read_VARS_LOOP_
	% read out parameters 
_Read_PARAMS_LOOP_
	% perform direct function calls
_DIR_FUN_CALL_LOOP_
	% evaluate the function values  
_Calc_EQUAS_AE_LOOP_
end

This code creates a Matlab function that reads all variable and parameter values, calls the functions of the model, and then calculates the equations of the model. Three subnodes are required (the block _DIR_FUN_CALL_LOOP_ is automatically included for the UDLS). Create the first subnode as “Variables”, block, the second one as “Parameters” block, and the third one as “EquationsAE” block. In the first block, add the following code:

	_VAR_NAME_ = X(_CTRVALUE_);

This code is again executed in a loop and writes the current entry of X into the variable name _VAR_NAME_. In the second block, add

	_PAR_NAME_ = PARAMS(_CTRVALUE_);

This does the same for the parameters. Finally, add the following code to the third block:

	Y(_CTRVALUE_) = _CALC_EXPR_;

This loop writes each equation into one entry of the variable Y.

Direct function calls

Now, we need to make sure that the functions within the model are called and implemented correctly. We start with the former by highlighting the block _DIR_FUN_CALL_LOOP_ and adding the following code:

	_FCT_OUT_VAR_NAME_ = _FUN_NAMESHORT_(_FUN_CALL_VARS_LOOP__SEP_SYMBOL__FUN_CALL_PARAMS_LOOP_);

This code specifies the output varnames of the functions within a loop. The functions are called with their shortnames with the call arguments and call parameters, both separated by the _SEP_SYMBOL_ (“,”).

Afterwards, we specify how the functions are implemented. Therefore, we highlight the block _DIR_FUN_IMPL_LOOP and add the code below:

function[_FCT_OUT_VAR_NAME_] = _FUN_NAMESHORT_(_FUN_IMPL_VARS_LOOP__SEP_SYMBOL__FUN_IMPL_PARS_LOOP_)
	_FCT_OUT_VAR_NAME_ = _FUN_EXPR_;
end

This code creates Matlab functions for all functions that are part of the model. The function returns the output variable name, is named according to its shortname (attention, here function call and function implementation must be consistent!), and receives the list of input variables and parameters, separated by the _SEP_SYMBOL_. Afterwards, the output variable is calculated using the function expression _FUN_EXPR_.

Display results

Finally, the solution for all iterated variables shall be displayed. Therefore, we add the block _Display_Results_ and add the following code:

function displayResults(SOL)

	% print variable values to display
_Disp_VARS_LOOP_
end

Within this function, we place another variable loop block with the following code:

	disp(['_VAR_NAME_ 	 ', num2str(SOL(_CTRVALUE_),'%.8f')]);

This code displays the variable name and the respective entry of the solution vector SOL with 8 decimal digits.

Now save the UDLS at a location of your preference.

Code generation

As already alluded to above, you can now go to the tab Overall Results, load a simulation, and click on Generate Full Code. This is good for debugging when you develop your UDLS. 

Once you are happy with the results, you can use the UDLS in the “Simulation” or “Optimization” section of MOSAICmodeling (here, we stick to simulations). In this example, we will load the simulation from the example Use of Functions III – Parameter Lists and Indices with ID 182676. Go to the Evaluation tab and then to the Generation tab. Within the Language Specification, select “User-defined” and open your saved UDLS. Note that you can now set the options you added in the tab Code Generation Properties. Once you are happy with the specificed options, click on Generate Code. The result is shown below. This code can now be copied to Matlab to solve the equation system.

% Matlab code to solve the nonlinear equation system below
% Code generated with MOSAICmodeling

function InitializeAndSolveEquationSystem()

	% load variable start values
	X_INIT(1) = 0.5;	% e0_L
	X_INIT(2) = 0.5;	% e0_V
	X_INIT(3) = 0.5;	% e0_x_c1
	X_INIT(4) = 0.5;	% e0_x_c2
	X_INIT(5) = 0.5;	% e0_y_c1
	X_INIT(6) = 0.5;	% e0_y_c2

	% load parameters
	PARAMS(1) = 1.0;	% e0_F
	PARAMS(2) = 354.6;	% e0_T
	PARAMS(3) = 0.5;	% e0_z_c1
	PARAMS(4) = 0.5;	% e0_z_c2
	PARAMS(5) = 101325.0;	% e0_P
	PARAMS(6) = 82.718;	% e0_A_c1
	PARAMS(7) = 73.649;	% e0_A_c2
	PARAMS(8) = -6904.5;	% e0_B_c1
	PARAMS(9) = -7258.2;	% e0_B_c2
	PARAMS(10) = -8.8622;	% e0_C_c1
	PARAMS(11) = -7.3037;	% e0_C_c2
	PARAMS(12) = 7.4664E-6;	% e0_D_c1
	PARAMS(13) = 4.1653E-6;	% e0_D_c2
	PARAMS(14) = 2.0;	% e0_E_c1
	PARAMS(15) = 2.0;	% e0_E_c2

	options = optimoptions('fsolve','MaxIter',1000,'MaxFunEvals',100000,'TolFun',1e-10,'Display','Iter');
	[SOL,RES,EXITFLAG] = fsolve (@( X_ITER ) modelEquations(X_ITER,PARAMS),X_INIT,options);

	displayResults(SOL);
	fprintf('\n')
	fprintf('Residuals: \n')
	disp(RES')
	fprintf('Terminated with exitflag %i \n',EXITFLAG)

end

function[Y] = modelEquations(X,PARAMS)

	% read out variables  
	e0_L = X(1);
	e0_V = X(2);
	e0_x_c1 = X(3);
	e0_x_c2 = X(4);
	e0_y_c1 = X(5);
	e0_y_c2 = X(6);

	% read out parameters 
	e0_F = PARAMS(1);
	e0_T = PARAMS(2);
	e0_z_c1 = PARAMS(3);
	e0_z_c2 = PARAMS(4);
	e0_P = PARAMS(5);
	e0_A_c1 = PARAMS(6);
	e0_A_c2 = PARAMS(7);
	e0_B_c1 = PARAMS(8);
	e0_B_c2 = PARAMS(9);
	e0_C_c1 = PARAMS(10);
	e0_C_c2 = PARAMS(11);
	e0_D_c1 = PARAMS(12);
	e0_D_c2 = PARAMS(13);
	e0_E_c1 = PARAMS(14);
	e0_E_c2 = PARAMS(15);

	% perform direct function calls
	e0_P_LV_c2 = fun_182674(e0_T,e0_A_c2,e0_B_c2,e0_C_c2,e0_D_c2,e0_E_c2);
	e0_P_LV_c1 = fun_182674(e0_T,e0_A_c1,e0_B_c1,e0_C_c1,e0_D_c1,e0_E_c1);

	% evaluate the function values  
	Y(1) = (e0_F)*(e0_z_c1)-((e0_L)*(e0_x_c1)+(e0_V)*(e0_y_c1));
	Y(2) = (e0_F)*(e0_z_c2)-((e0_L)*(e0_x_c2)+(e0_V)*(e0_y_c2));
	Y(3) = (e0_x_c1)*(e0_P_LV_c1)-((e0_y_c1)*(e0_P));
	Y(4) = (e0_x_c2)*(e0_P_LV_c2)-((e0_y_c2)*(e0_P));
	Y(5) = 1.0-((e0_x_c1+e0_x_c2));
	Y(6) = 1.0-((e0_y_c1+e0_y_c2));

end

function[std_P_LV] = fun_182674(std_T,std_A,std_B,std_C,std_D,std_E)
	std_P_LV = exp(std_A+(std_B)/(std_T)+(std_C)*(log(std_T))+(std_D)*(((std_T))^(std_E)));
end

function displayResults(SOL)

	% print variable values to display
	disp(['e0_L 	 ', num2str(SOL(1),'%.8f')]);
	disp(['e0_V 	 ', num2str(SOL(2),'%.8f')]);
	disp(['e0_x_c1 	 ', num2str(SOL(3),'%.8f')]);
	disp(['e0_x_c2 	 ', num2str(SOL(4),'%.8f')]);
	disp(['e0_y_c1 	 ', num2str(SOL(5),'%.8f')]);
	disp(['e0_y_c2 	 ', num2str(SOL(6),'%.8f')]);

end