mp.htm

optimization toolbox

      </table>
      <p>Now, instead of setting up the problem for the whole prediction 
		horizon, we only set it up for one step, solve the problem 
		parametrically, take one step back, and perform a standard dynamic 
		programming value iteration.</p>
      <table cellPadding="10" width="100%" id="table13">
        <tr>
          <td class="xmpcode">
          <pre>J{N} = 0;

for k = N-1:-1:1    
    
    % Feasible region
    F =     set(-1 &lt; u{k}     &lt; 1);
    F = F + set(-1 &lt; C*x{k}   &lt; 1);
    F = F + set(-5 &lt; x{k}     &lt; 5);
    F = F + set(-1 &lt; C*x{k+1} &lt; 1);
    F = F + set(-5 &lt; x{k+1}   &lt; 5);   </pre>
			<pre>    % Dynamics
    F = F + set(x{k+1} == A*x{k}+B*u{k});</pre>
			<pre>    % Cost in value iteration
    obj = norm(x{k},1) + norm(u{k},1) + J{k+1}

    % Solve one-step problem    
    [sol{k},diagnost{k},Uz{k},J{k},Optimizer{k}] = solvemp(F,obj,[],x{k},u{k});
end</pre>
          </td>
        </tr>
      </table>
      <p>Notice the minor changes needed compared to the one shot solution. 
		Important to understand is that the value function at step <b>k</b> will 
		be a function in <b>x{k}</b>, hence when it is used at <b>k-1</b>, it 
		will be a function penalizing the predicted state. Note that YALMIP 
		automatically keep track of convexity of the value function. Hence, for 
		this example, no binary variables are introduced along the solution 
		process (this is why we safely can 
		omit the use of <a href="reference.htm#bounds">bounds</a>).</p>
		<p>To study the development of the value function, we can easily plot 
		them. </p>
      <table cellPadding="10" width="100%" id="table14">
        <tr>
          <td class="xmpcode">
          <pre>for k = N-1:-1:1
 plot(J{k});hold on
end</pre>
          </td>
        </tr>
      </table>
      <p>The optimal control input can also be plotted.</p>
      <table cellPadding="10" width="100%" id="table15">
        <tr>
          <td class="xmpcode">
          <pre>clf;plot(Optimizer{1})</pre>
          </td>
        </tr>
      </table>
      <p>Of course, you can do the same thing by working directly with the 
		numerical MPT 
		data.</p>
      <table cellPadding="10" width="100%" id="table23">
        <tr>
          <td class="xmpcode">
          <pre>for k = N-1:-1:1
 mpt_plotPWA(sol{k}{1}.Pn,sol{k}{1}.Bi,sol{k}{1}.Ci);hold on
end</pre>
          </td>
        </tr>
      </table>
      <p>Why not solve the problem with a worst-case L<sub><font face="Arial">&#8734;</font></sub> cost! 
		Notice&nbsp;the non-additive value iteration! Although this might look 
		complicated using a mix of maximums, norms and piecewise value-functions 
		, the formulation fits perfectly into the <a href="extoperators.htm">nonlinear operator</a> 
		framework in YALMIP.</p>
      <table cellPadding="10" width="100%" id="table24">
        <tr>
          <td class="xmpcode">
          <pre>J{N} = 0;

for k = N-1:-1:1    
    
    % Feasible region
    F =     set(-1 &lt; u{k}     &lt; 1);
    F = F + set(-1 &lt; C*x{k}   &lt; 1);
    F = F + set(-5 &lt; x{k}     &lt; 5);
    F = F + set(-1 &lt; C*x{k+1} &lt; 1);
    F = F + set(-5 &lt; x{k+1}   &lt; 5);   </pre>
			<pre>    % Dynamics
    F = F + set(x{k+1} == A*x{k}+B*u{k});</pre>
			<pre>    % Cost in value iteration
    obj = max(norm([x{k};u{k}],inf),J{k+1})

    % Solve one-step problem    
    [sol{k},diagnost{k},Uz{k},J{k},Optimizer{k}] = solvemp(F,obj,[],x{k},u{k});
end</pre>
          </td>
        </tr>
      </table>
      <h3>Dynamic programming with PWA systems</h3>
		<p>As mentioned above, the difference between the value function from a 
		parametric LP and a parametric LP with binary variables is that 
		convexity of the value function no longer is guaranteed. In practice, 
		this means that (additional) binary variables are required when the the 
		value function is used in an optimization problem. This is done 
		automatically in YALMIP through the <a href="extoperators.htm">nonlinear 
		operator</a> framework, but it is extremely important to realize that 
		the value function from a binary problem is much more complex than the 
		one resulting from a standard parametric problem. Nevertheless, working 
		with the objects is easy, as the following example illustrates. In this 
		case, we solve the dynamic programming problem for our PWA system</p>
      <table cellPadding="10" width="100%" id="table20">
        <tr>
          <td class="xmpcode">
          <pre>yalmip('clear')
clear all</pre>
			<pre>% Model data
A = [2 -1;1 0];
B1 = [1;0];
B2 = [pi;0]; % Larger gain for negative first state
C = [0.5 0.5];</pre>
			<pre>nx = 2; % Number of states
nu = 1; % Number of inputs

% Prediction horizon
N = 4;</pre>
          <pre>% States x(k), ..., x(k+N)
x = sdpvar(repmat(nx,1,N),repmat(1,1,N));</pre>
			<pre>% Inputs u(k), ..., u(k+N) (last one not used)
u = sdpvar(repmat(nu,1,N),repmat(1,1,N));</pre>
			<pre>% Binary for PWA selection
d = binvar(repmat(2,1,N),repmat(1,1,N));</pre>
          </td>
        </tr>
      </table>
      <p>Just as in the LTI case, we set up the problem for the one step case, 
		and use the value function from the previous iteration.</p>
      <table cellPadding="10" width="100%" id="table21">
        <tr>
          <td class="xmpcode">
          <pre>J{N} = 0;

for k = N-1:-1:1    

    % Strenghten big-M (improves numerics)
    bounds(x{k},-5,5);
    bounds(u{k},-1,1);
    bounds(x{k+1},-5,5);
    
    % Feasible region
    F =     set(-1 &lt; u{k}     &lt; 1);
    F = F + set(-1 &lt; C*x{k}   &lt; 1);
    F = F + set(-5 &lt; x{k}     &lt; 5);
    F = F + set(-1 &lt; C*x{k+1} &lt; 1);
    F = F + set(-5 &lt; x{k+1}   &lt; 5);  </pre>
			<pre>    % PWA Dynamics
    F = F + set(implies(d{k}(1),x{k+1} == (A*x{k}+B1*u{k})));
    F = F + set(implies(d{k}(2),x{k+1} == (A*x{k}+B2*u{k})));
    F = F + set(implies(d{k}(1),x{k}(1) &gt; 0));
    F = F + set(implies(d{k}(2),x{k}(1) &lt; 0));
    F = F + set(sum(d{k}) == 1);</pre>
			<pre>    % Cost in value iteration
    obj = norm(x{k},1) + norm(u{k},1) + J{k+1}</pre>
			<pre>    % Solve one-step problem
    [sol{k},diagnost{k},Uz{k},J{k},Optimizer{k}] = solvemp(F,obj,[],x{k},u{k});
end</pre>
          </td>
        </tr>
      </table>
      <p>Plotting the final value function clearly reveals the nonconvexity.</p>
      <table cellPadding="10" width="100%" id="table22">
        <tr>
          <td class="xmpcode">
          <pre>clf;plot(J{1})</pre>
          </td>
        </tr>
      </table>
      <h3>Behind the scenes and advanced use</h3>
		<p>The first thing that might be a bit unusual to the advanced user is 
		the piecewise functions that YALMIP returns in the fourth and fifth 
		output from <code>solvemp</code> . In principle, they are specialized
		<a href="reference.htm#pwf">pwf</a> objects. To create a PWA value 
		function after solving a multi-parametric LP, the following 
		command is used.</p>
      <table cellPadding="10" width="100%" id="table26">
        <tr>
          <td class="xmpcode">
          <pre>Valuefunction = pwf(sol,x,'convex')</pre>
          </td>
        </tr>
      </table>
      <p>The <a href="reference.htm#pwf">pwf</a> command will recognize the MPT 
		solution structure and create a PWA function based on the 
		fields <b>Pn</b>, <b>Bi</b> and <b>Ci</b>. The dedicated
		<a href="extoperators.htm">nonlinear operator</a> is implemented in the file <code>pwa_yalmip.m</code> 
		. The <a href="extoperators.htm">nonlinear operator</a> will exploit the 
		fact that the PWA function is convex and implements an efficient epi-graph 
		representation. In case the PWA function is used in a non-convex fashion 
		(i.e. the automatic <a href="extoperators.htm#propagation">convexity 
		propagation</a> fails), a
		<a href="extoperators.htm#milp">MILP implementation</a> is also 
		available.<p>If the field <b>Ai</b> is non-empty (solution obtained from 
		a multi-parametric QP), a corresponding PWQ function is created (<code>pwq_yalmip.m</code> 
		). <p>To create a PWA function 
		representing the optimizer, two things have to be changed. To begin 
		with, YALMIP searches for the <b>Bi</b> and <b>Ci</b> fields, but since 
		we want to create a PWA function based on <b>Fi</b> and <b>Gi</b> 
		fields, the field names have to be changed. Secondly, the piecewise&nbsp; 
		optimizer is typically not convex, so a general PWA function is created 
		instead (requiring 1 binary variable per region if the variable later 
		actually is used in an optimization problem.)<table cellPadding="10" width="100%" id="table27">
        <tr>
          <td class="xmpcode">
          <pre>sol.Ai = cell(1,length(sol.Ai));
sol.Bi = sol.Fi;
sol.Ci = sol.Gi;
Optimizer = pwf(sol,x,'general')</pre>
          </td>
        </tr>
      </table>
      <p>A third important case is when the solution structure returned from <code>solvemp</code> is 
		a cell with several <a href="solvers.htm#mpt">MPT</a> structures. This means that a 
		multi-parametric problem with binary variables was 
		solved, and the different cells represent overlapping solutions. One way 
		to get rid of the overlapping solutions is to use the MPT command 
		<code>mpt_removeOverlaps.m</code> and create a PWA function based on the result. Since the 
		resulting PWA function typically is non-convex, 
		we must create a general function.<table cellPadding="10" width="100%" id="table28">
        <tr>
          <td class="xmpcode">
          <pre>sol_single = mpt_removeOverlaps(sol);
Valuefunction = pwf(sol_single,x,'general')</pre>
          </td>
        </tr>
      </table>
		<p>This is only recommended if you just intend to plot or investigate 
		the value function. Typically, if you want to continue computing with the 
		result, i.e. use it in another optimization problem, as in the dynamic 
		programming examples above, it is recommended to keep the overlapping 
		solutions. To model a general PWA function, 1 binary variable per region 
		is needed. However, by using a dedicated <a href="extoperators.htm">nonlinear operator</a> 
		for overlapping convex PWA functions, only one binary per PWA function 
		is needed.</p>
		<table cellPadding="10" width="100%" id="table29">
        <tr>
          <td class="xmpcode">
          <pre>Valuefunction = pwf(sol,x,'overlapping')</pre>
          </td>
        </tr>
      </table>
		<p>As we have seen in the examples above, the PWA and PWQ functions can be 
		plotted. This is currently nothing but a simple hack. Normally, when you 
		apply the plot command, the corresponding double values are plotted. 
		However, if the input is a simple scalar PWA or PWQ variable, the 
		underlying MPT structures will be extracted and the plot commands in MPT 
		will be called.</td>
    </tr>
  </table>
	<p>&nbsp;</div>

</body>

</html>