Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I write supplemental explanation to render my aim of optimization.

First Consider discrete time linear system below:

\begin{equation} G(z)=\frac{1-e^{-0.3T}}{z^{100}(z-e^{-0.3T})} \end{equation}

where $z = e^{j\omega}$ and $\omega$ is frequency, $\omega \in [-\pi,\pi)$, and $T$ is sampling period, assume that $T = 1$. It's obvious that the absolute value of $z$ is equal to $1$ for all $\omega$ between $-\pi$ and $\pi$.

When one is talking about maximum absolute value of $G(z)$, in fact he/she consider below function :

\begin{equation} \max_{\omega} |G(e^{j\omega})|=|\frac{1-e^{-0.3T}}{e^{j\omega}-e^{-0.3T}}| . \end{equation}

note that $|z^{100}| = |z|^{100} = 1.$

On the other hand consider another discrete time linear system :

\begin{equation} C(z) = \frac{(x+\frac{T}{2}y)z+(\frac{T}{2}y-x)}{z-1} \end{equation}

Consider these two linear system connected serially.Then the closed loop system transfer function (with unity feedback) is:

\begin{equation} f(x,y,z) = \frac{G(z)C(z)}{1+G(z)C(z)}\end{equation}


\begin{equation} f(x,y,e^{j\omega}) = \frac{G(e^{j\omega})C(e^{j\omega})}{1+G(e^{j\omega})C(e^{j\omega})}\end{equation}

My aim is to minimize the maximum absolute value of $f$ (respect to $\omega$) over $x $ and $ y$ as mathematically described bellow:

\begin{equation} \min_{x, y} \max_{z \; or \; \omega} |f(x,y,z)|\end{equation}

This is equal to minimizing the $L_\infty$ norm of $f$ (respect to $\omega$) over $x$ and $y$ (right?) In fact:

\begin{equation}\min_{x, y} \max_{\omega} |f(x,y,\omega)| \equiv \min_{x, y} ||f(x,y,\omega)||_{\infty , \;\omega}\end{equation}

How can I do the optimization with Mathematica?

I try the code below for defining the function and maximizing $f$ respect to $z$:

T = 0.01
z = Exp[I w]
f = (((x + (T/2) y) z + ((T/2) y - x)) (1 - Exp[-0.3 T]))/(
      (z^100) (z - Exp[-0.3 T]) (z - 1) + ((x + (T/2) y) z + ((T/2) y - x)) (1 - Exp[-0.3 T]))
maxp = First[ NMaximize[{f[x, y, w], 0 <= w <= 10}, w, Method -> DifferentialEvolution]]

But when I run code, the following error appears:

NMaximize::nnum: The function value 
-((0.0029955 (-x + (0.89724 + 0.441543 I) (x + 0.005 y) + 0.005 y))/((0.120832 -0.165872 I)
 + 0.0029955 ( -x + (0.89724 + <<20>> I) (<<1>>)+0.005 y)))[x,y,0.457318] 
is not a number at {w} = {0.457318}. >>

I try different intervals for $w$ but nothing changed.

What's wrong with it?

share|improve this question
This can be done as a bilevel optimization. See this MathGroup post or this Matlab Usenet group post or the "fewnomial" example here – Daniel Lichtblau Jan 22 '14 at 23:19
@DanielLichtblau Thanks. But I have problem with doing that. I added more details in question. – Zia Jan 23 '14 at 10:47
That's very confused code. First, you might want f[x_,y_,z_]:=... Then there is the bit about defining z=Exp[I*w] and later invoking f[x,y,w]. Very confusing as to what specifically you want for that third argument. Finally there is absolutely no indication of what you do to minimize over {x,y} , and that is precisely the sort of thing covered at the links I posted. – Daniel Lichtblau Jan 23 '14 at 15:51
@DanielLichtblau I'm new with mathematica. Sorry for these silly mistakes. I don't know what do you mean by doing " [x_,p_Real,n_]". Could you explain more? – Zia Jan 23 '14 at 16:10
A pattern such as _Real in e.g f[x_,p_Real,n_] means the function f[] will not evaluate UNLESS that argument is an explicit real number, e.g. f[x,5.3,2]. This is useful for preventing a calling function like` NMinimize` from treating f[x,p,n] as symnolic, attempting preprocessing to find say a gradient, that sort of thing. – Daniel Lichtblau Jan 23 '14 at 16:46
up vote 10 down vote accepted

This is along the right lines, I think. But I had to change the minimization to use a log, in order for the scaling not to be too much of a problem for FindMinimum. I also removed the complex exponential definition of z. You might want to exponentiate the result. I also used an Abs to make the logarithm real valued, so that too could have consequences.

t = 0.01;
f[x_, y_, z_] := 
 Log[Abs[(((x + (t/2) y) z + ((t/2) y - x)) (1 - 
        Exp[-0.3 t]))/((z^100) (z - Exp[-0.3 t]) (z - 
         1) + ((x + (t/2) y) z + ((t/2) y - x)) (1 - Exp[-0.3 t]))]]

fzmin[z_?NumberQ] := Module[{x, y, res},
  res = FindMinimum[f[x, y, z], {x, -10, 10}, {y, -10, 10}];
  If[! ListQ[res], -1000000., res[[1]]]]

max = 
 NMaximize[{fzmin[z], 0 <= z <= 10}, {z, 0, 10}, 
  Method -> "DifferentialEvolution"]

(* Out[317]= {-1.11022302463*10^-16, {z -> 0.416231991971}} *)

--- edit ---

Based on updates to the question, the code below should be appropriate. I remark that it is not hugely modified from the prior attempt.

T = 0.01;
f[x_, y_, w_] := (((x + (T/2) y) Exp[I w] + ((T/2) y - x)) (1 - 
       Exp[-0.3 T]))/((Exp[I w]^100) (Exp[I w] - 
        Exp[-0.3 T]) (Exp[I w] - 1) + 
        ((x + (T/2) y) Exp[I w] + ((T/2) y - x)) (1 - Exp[-0.3 T]));

g[x_?NumberQ, y_?NumberQ] := 
 NMaximize[{Evaluate[Abs[f[x, y, w]]], 0 <= w <= 10}, {w, 0, 10}][[1]]

NMinimize[g[x, y], {x, y}]

{0.777743387928, {x -> -6.56647441199, y -> -2.87876723563}}

--- end edit ---

share|improve this answer
But i doubt on the way you propose. why you first minimize $f$ then maximized the result? This is not compatible with question. I added more details to question. $z$ should be equal to $e^{j\omega}$. Thanks for your answer. – Zia Jan 24 '14 at 9:05
Thanks. I run the above code but process time is too long! One question: Why you use [[1]] in code? How does it help? – Zia Jan 24 '14 at 19:03
(1) I would not expect a bilevel optimization to be fast. But you might be able to improve the speed by using FindMaximum and FindMinimum in place of NMaximize and NMinimize. (2) It extracts the first component of what NMaximize returns. That's what you want. The max value. Not the replacement rules that give the arg max. (3) You'll need to learn a bit more about Mathematica in order to use it effectively. I'd recommend experimenting some. For example, you could remove that [[1]] and run e.g. g[1.2,3.3] to see what it then does differently. – Daniel Lichtblau Jan 24 '14 at 21:15
Thanks for everything. I hope to work more with mathematica and ask better questions in the future. – Zia Jan 25 '14 at 4:38

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.