Use of User-defined Language Specifications II – Dynamic Optimization via Single Shooting

Goal

In this example, we will set up a UDLS for sequential optimization in Matlab. The UDLS is available in MOSAICmodeling with ID 185804.

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. Deactivate Direct Functions Supported and Case Sensitive in the Language Settings. The former will prohibit the inclusion of direct function calls, the latter ensures a difference between capital and small letters. Note that we exclude functions only for the sake of simplicity and their addition to the UDLS would be fairly simple.
  2. Activate Force Floating Point in the Number settings
  3. Select “DAE” as Equation System Type
  4. Select “Optimization” and “Sequential” in the box Opt or Nopt

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.
  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 this example and can therefore be neglected.

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

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 to add options for the user export. Which options exactly you would like to define, is up to the user. In the UDLS, there is an extensive list of options for fmincon (the optimizer), DAE solvers (for integrating the dynamic model), and for fsolve (to solve the algebraic sub-system). If you want to see the complete list, take a look at the UDLS with the ID given above. Here, we just want to highlight a couple of options as they are particularly relevant for the example Use of Optimizations II – Dynamic Optimization via Single Shooting, in which this UDLS is applied.

  1. Add the the Property OdeSolver with Handle _ODESOLVER_, Type STRING, and Default ode15
  2. Add the the Property Algorithm with Handle _ALGORITHM_, Type STRING, and Default ‘sqp’
  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 MaxIterations with Handle _MAXITERATIONS_, Type INTEGER, and Default 100000
  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’, ‘notify-detailed’, ‘iter-detailed’, and ‘final-detailed’. This allows the user to select a preferred choice within MOSAICmodeling when a simulation shall be exported

Code elements

This tab is fairly complicated as it contains the structure for the generated code. 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_

_Optimization_layer_

_Simulation_layer_

_Model_Equations_

_DIR_FUN_IMPL_LOOP_

_Display_Results_

These blocks will execute the following tasks:

  • _Help_Header_: generation of documentation code at the beginning
  • _Optimization_layer_: variable initialization and call of fmincon
  • _Simulation_layer_: integration of DAE system and calculation of objective function
  • _Model_Equations_: model equations of DAE system, algebraic sub-system, and nonlinear constraints for optimization
  • _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 optimize the differential-algebraic 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.

Optimization layer

In the optimization layer, the variables are initialized, the decisions and their bounds are set, the solver options are specified, and the optimizer is called according to the code below:

function SolveDynamicOptimizationProblem()
	clear
	close all
	clc

	% declare interval of independent variable
_Set_DIFF_VAR_
	% initial conditions for explicit state variables
_Set_Init_STATE_VARS_EXPL_LOOP_
	% initial conditions for implicit state variables
_Set_STATE_VARS_IMPL_LOOP_
	% objective function
_Set_OBJECTIVE_LOOP_
	% set parameters
_Set_PARAMS_LOOP_
	% define decision variables and their bounds
_Set_DECISIONS_
	M = [eye(_DIM_DIFFEQUAS_,_DIM_DIFFEQUAS_),zeros(_DIM_DIFFEQUAS_,_DIM_ALGEQUAS_); zeros(_DIM_ALGEQUAS_, _DIM_ALGEQUAS_+_DIM_DIFFEQUAS_)];

	% specify options	
	OPToptions = optimoptions('fmincon',...
		'Algorithm',_ALGORITHM_,...
		'CheckGradients',_CHECKGRADIENTS_,...
		'ConstraintTolerance',_CONSTRAINTTOLERANCE_,...
		'Diagnostics',_DIAGNOSTICS_,...
		'DiffMaxChange',_DIFFMAXCHANGE_,...
		'DiffMinChange',_DIFFMINCHANGE_,...
		'Display',_DISPLAY_,...
		'FiniteDifferenceStepSize',_FINITEDIFFERENCESTEPSIZE_,...
		'FiniteDifferenceType',_FINITEDIFFERENCETYPE_,...
		'FunValCheck',_FUNVALCHECK_,...
		'MaxFunctionEvaluations',_MAXFUNCTIONEVALUATIONS_,...
		'MaxIterations',_MAXITERATIONS_,...
		'OptimalityTolerance',_OPTIMALITYTOLERANCE_,...
		'OutputFcn',_OUTPUTFCN_,...
		'PlotFcn',_PLOTFCN_,...
		'SpecifyConstraintGradient',_SPECIFYCONSTRAINTGRADIENT_,...
		'SpecifyObjectiveGradient',_SPECIFYOBJECTIVEGRADIENT_,...
		'StepTolerance',_STEPTOLERANCE_,...
		'TypicalX',_TYPICALX_,...
		'UseParallel',_USEPARALLEL_);

	DAEoptions = odeset('AbsTol',_DAEABSTOL_,...
		'BDF',_DAEBDF_,...
		'Events',_DAEEVENTS_,...
		'InitialStep',_DAEINITIALSTEP_,...
		'Jacobian',_DAEJACOBIAN_,...
		'JConstant',_DAEJCONSTANT_,...
		'JPattern',_DAEJPATTERN_,...
		'Mass',M,...
		'MassSingular',_DAEMASSSINGULAR_,...
		'MaxOrder',_DAEMAXORDER_,...
		'MaxStep',_DAEMAXSTEP_,...
		'MStateDependence',_DAEMSTATEDEPENDENCE_,...
		'MvPattern',_DAEMVPATTERN_,...
		'NonNegative',_DAENONNEGATIVE_,...
		'NormControl',_DAENORMCONTROL_,...
		'OutputFcn',_DAEOUTPUTFCN_,...
		'OutputSel',_DAEOUTPUTSEL_,...
		'Refine',_DAEREFINE_,...
		'RelTol',_DAERELTOL_,...
		'Stats',_DAESTATS_,...
		'Vectorized',_DAEVECTORIZED_);

	NLEoptions = optimoptions('fsolve',...
		'Display',_NLEDISPLAY_,...
		'FunctionTolerance',_NLEFUNCTIONTOLERANCE_,...
		'MaxFunctionEvaluations',_NLEMAXFUNCTIONEVALUATIONS_,...
		'MaxIterations',_NLEMAXITERATIONS_,...
		'OptimalityTolerance',_NLEOPTIMALITYTOLERANCE_);

	% solve sequential optimization problem
	[DEC,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = fmincon(@(DEC) objectiveFunction(DEC,X_INTERVAL,Y_INIT,PARAMS,NLEoptions,DAEoptions),DEC_INIT,[],[],[],[],LB,UB,@(DEC) NONLCON(DEC,X_INTERVAL,Y_INIT,PARAMS,NLEoptions,DAEoptions),OPToptions);

	% plot optimized profiles
	[X,Y] = ode15s(@(X,Y)DAESystem(X,Y,PARAMS,DEC),X_INTERVAL,Y_INIT',DAEoptions);
	displayResults(X,Y);

end

This code creates a function in Matlab that begins with a looping setting the time horizon of the integration (which is in this case equal to the optimization horizon), set the initial guesses for both explicit and implicit variables. Finally, the objective is – in this specific case – added to the variable vector of the DAE system. Then, the parameter values are set within a loop. Finally, the decisions and their bounds are set.

We then proceed with the specification of the mass matrix for the DAE solver and the various options for optimizer, DAE solver, and NLE solver. Then, we call fmincon to solve the optimization problem. Finally, the DAE system is solved again with the decisions obtained via optimization and the profiles of the variables are plotted.

As you can see, the block requires the loops _Set_DIFF_VAR_ (DifferentialVariable), _Set_STATE_VARS_EXPL_LOOP_ (StateVariablesExplicit), _Set_STATE_VARS_IMPL_LOOP_ (StateVariablesImplicit), _Set_OBJECTIVE_LOOP_ (ObjectiveFunction), _Set_DECISIONS_ (GenericVariables), and _Set_PARAMS_LOOP_ (Parameters). These loops are now added as subnodes of the node _Optimization_layer_. For most of these nodes, the code to be entered is straightforward, e.g., for the parameter block:

	Y_INIT(_CTRVALUE_) = _INIT_VALUE_;  	% _VAR_NAME_

For the GenericVariables block, a configuration is required. Therefore, right-click on the node _Set_Decisions_ and set OPTIMIZATION_VAR to True. Then add the code below to the block:

	DEC_INIT(_CTRVALUE_) = _VAR_VALUE_;	% _VAR_NAME_
	LB(_CTRVALUE_) = _LOWER_BOUND_;
	UB(_CTRVALUE_) = _UPPER_BOUND_;

This ensures that only the degrees of freedom of the optimization problem are initialized with value and bounds.

Simulation layer

After you have added code to all blocks described above, create a new sub-node of the Full Code and name it _Simulation_layer. Enter the code below into this block:

function OBJ = objectiveFunction(DEC,X_INTERVAL,Y_INIT,PARAMS,NLEoptions,DAEoptions)

	% if true, the algebraic sub-system of the DAE is solved
	if _SOLVEALGEBRAICSUBSYSTEM_
		Y_DIFF_INIT = Y_INIT(1:_DIM_DIFFEQUAS_);
		Y_ALG_INIT = Y_INIT(1+_DIM_DIFFEQUAS_:end);
		Y_ALG_INIT = fsolve (@( Y_ALG ) AESubSystem(Y_ALG,Y_DIFF_INIT,PARAMS,DEC),Y_ALG_INIT,NLEoptions);
		Y_INIT = [Y_DIFF_INIT, Y_ALG_INIT];
	end

	[X,Y] = _ODESOLVER_(@(X,Y)DAESystem(X,Y,PARAMS,DEC),X_INTERVAL,Y_INIT',DAEoptions);

_OBJECTIVE_LOOP_

end

This code optionally solves the algebraic sub-system to initialize the DAE and then calls the DAE solver specified by the user to integrate the system. Within the code, the _OBJECTIVE_LOOP_ is a sub-node of type ObjectiveFunction with the following code:

	OBJ = Y(end,end);

This code applies if we are interested in the value of the objective at the end. In our application example, the objective is equal to the negative integral of the molar concentration of a specific component. As the integral is formulated as differential equation, we simply need its value at the end of the time horizon.

Model equations

The next block writes the DAE system, the algebraic sub-system, and the nonlinear (in)equality constraints:

 function DYDX = DAESystem(X,Y,PARAMS,DEC)

	% read out differential variable
_Read_DIFF_VAR_
	% read out parameters 
_Read_PARAMS_LOOP_
	% read out explicit state variables  
_Read_STATE_VARS_EXPL_LOOP_
	% read out implicit state variables
_Read_STATE_VARS_IMPL_LOOP_
_Read_OBJECTIVE_LOOP_
_Read_DECISIONS_
_DIR_FUN_CALL_LOOP_
	% evaluate the explicit differential equations 
_Calc_EQUAS_ODE_LOOP_
	% evaluate the implicit algebraic equations
_Calc_EQUAS_AE_LOOP_
	DYDX=DYDX';

end

function DYDX = AESubSystem(Y_ALG,Y_DIFF,PARAMS,DEC)

	Y = [Y_DIFF, Y_ALG];

	% read out parameters 
_Read_PARAMS_LOOP_
	% read out explicit state variables  
_Read_STATE_VARS_EXPL_LOOP_
	% read out implicit state variables
_Read_STATE_VARS_IMPL_LOOP_
_Read_OBJECTIVE_LOOP_
_Read_DECISIONS_
_DIR_FUN_CALL_LOOP_
	% evaluate the explicit differential equations 
_Calc_EQUAS_ODE_LOOP_
	% evaluate the implicit algebraic equations
_Calc_EQUAS_AE_LOOP_
	DYDX=DYDX(_DIM_DIFFEQUAS_+1:end)';

end

function [c, ceq] = NONLCON(DEC,X_INTERVAL,Y_INIT,PARAMS,NLEoptions,DAEoptions)

	% if true, the algebraic sub-system of the DAE is solved
	if _SOLVEALGEBRAICSUBSYSTEM_
		Y_DIFF_INIT = Y_INIT(1:_DIM_DIFFEQUAS_);
		Y_ALG_INIT = Y_INIT(1+_DIM_DIFFEQUAS_:end);
		Y_ALG_INIT = fsolve (@( Y_ALG ) AESubSystem(Y_ALG,Y_DIFF_INIT,PARAMS,DEC),Y_ALG_INIT,NLEoptions);
		Y_INIT = [Y_DIFF_INIT, Y_ALG_INIT];
	end

	[X,Y] = _ODESOLVER_(@(X,Y)DAESystem(X,Y,PARAMS,DEC),X_INTERVAL,Y_INIT',DAEoptions);

	% inequality constraints for lower bound
_Calc_Bounds_STATE_VARS_EXPL_LOOP_
_Calc_Bounds_STATE_VARS_IMPL_LOOP_
	c = [c_lb, c_ub];
	ceq = [];

end

Again, there are a couple of sub-nodes to read the specific variable and parameter names from the generic variable and parameter vectors Y and PARAMS, to read the decision variables DEC, etc. Details can be found in the UDLS with the ID given at the top of this page.

Display results

Finally, we want to be able to plot the dynamic profiles of the variables. For this purpose, we define the relatively simple plot function below. More advanced would be to also plot the decision variables over time, etc.

function displayResults(X,Y)
	% decide for a plot type: 
	%   0 	-> Plot the variables into individual figures 
	%   1 	-> Plot into sub figures 
	%   2 	-> Plot all selected into one figure 
	%   other -> Do not plot 
	plotType = _PLOTTYPE_; 

	% set a line width: 
	linewidth = 1.5; 

	% define which dependent variables should be plotted
	%   1 	-> Plot.
	%   other 	-> Do not plot.
_Plot_Control_State_Vars_

	% decide wether to normalize the y axis
	%   1 	-> Normalized
	%   other 	-> Individual maximum scale
	axisControl = 1;

	%====================================================

	% labels of the dependent variables
_Label_State_Vars_
_Label_Diff_Var_

	% plot the variables 
	figureIndex=1; 
	xMinVal = min(X); 
	xMaxVal = max(X); 
	yMinVal = min(min(Y)); 
	yMaxVal = max(max(Y)); 
	if (plotType==0) 
		% create a plot for each state variable individually 
		for i=1:_DIM_VARS_
			if (plotControl(i)==1)
				fig=figure(i);
				set(fig,'defaultTextInterpreter','none')
				plot(X,Y(:,i),'LineWidth',linewidth)
				title('Solution of the DAE equation system');
				xlabel(xAxisLabel);
				ylabel(yAxisLabels{i});
				if (axisControl==1)
					axis([xMinVal xMaxVal min(Y(:,i)) max(Y(:,i))]);
				end
				legend(yAxisLabels(i,:));
				figureIndex=figureIndex+1;
			end
		end
	elseif (plotType==1)
		% use a subplot environment
		firstTime = 1;
		numberOfPlots = sum(plotControl);
		fig=figure(1);
		set(fig,'defaultTextInterpreter','none')
		for i=1:_DIM_VARS_
			if (plotControl(i)==1)
				subplot(numberOfPlots,1,figureIndex);
				plot(X,Y(:,i),'LineWidth',linewidth)
				ylabel(yAxisLabels{i});
				legend(yAxisLabels{i});
				if (axisControl==1)
					axis([xMinVal xMaxVal min(Y(:,i)) max(Y(:,i))]);
				end
				figureIndex=figureIndex+1;
				if (firstTime) 
					title('Solution of the DAE equation system');
					firstTime = 0;
				end 
			end 
		end
		xlabel(xAxisLabel);
	elseif (plotType==2)
		% plot in one figure
		colors = [
			'r'	% -> Red
			'g' % -> Green
			'b' % -> Blue
			'c' % -> Cyan
			'm' % -> Magenta
			'y' % -> Yellow
			'k' % -> Black
			];
		colorCtr=1;
		maxColors=7;
		fig=figure(1);
		set(fig,'defaultTextInterpreter','none')
		hold on;
		for i=1:_DIM_VARS_
			if (plotControl(i)==1)
				linespec = colors(colorCtr);
				if (colorCtr<=maxColors)
					colorCtr = colorCtr+1;
				end
				plot(X,Y(:,i),linespec,'LineWidth',linewidth);
			end
		end
		title('Solution of the DAE equation system');
		xlabel(xAxisLabel);
		ylabel(yAxisLabels{i});
		legend(yAxisLabels{:});
		hold off;
	end


end

Once you have tested the UDLS, save it so that is available in the future.