Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
How do I set x0 in quadprog to get better solution?

Subject: How do I set x0 in quadprog to get better solution?

From: William

Date: 22 Mar, 2006 04:56:01

Message: 1 of 30

Hi everyone! I've recently used 'quadprog' function to solve an
optimization of a 300X300 Hessian matrix:
X=quadprog(H,f,[],[],Aeq,beq,LB,UB,X0), and Matlab 7 tends to switch to
a medium-scale method because of the restriction conditions. Yet the
process is slow and get a bad solution with exceeding the maximum
number of iterations. And Matlab 7 gives advices:
  'Maximum number of iterations exceeded; increase options.MaxIter.To
continue solving the problem with the current solution as the starting
point, set x0 = x before calling quadprog.'

  Considering the large Hessian matrix, I am unwilling to increase
options.MaxIter. But I am quite interested in whether I can get a
better solution by 'set x0 = x before calling quadprog' or not, and
how?

  Thank you very much!

Subject: How do I set x0 in quadprog to get better solution?

From: John D'Errico

Date: 22 Mar, 2006 08:21:05

Message: 2 of 30

William wrote:
>
>
> Hi everyone! I've recently used 'quadprog' function to solve an
> optimization of a 300X300 Hessian matrix:
> X=quadprog(H,f,[],[],Aeq,beq,LB,UB,X0), and Matlab 7 tends to
> switch to
> a medium-scale method because of the restriction conditions. Yet
> the
> process is slow and get a bad solution with exceeding the maximum
> number of iterations. And Matlab 7 gives advices:
> 'Maximum number of iterations exceeded; increase
> options.MaxIter.To
> continue solving the problem with the current solution as the
> starting
> point, set x0 = x before calling quadprog.'
>
> Considering the large Hessian matrix, I am unwilling to increase
> options.MaxIter. But I am quite interested in whether I can get a
> better solution by 'set x0 = x before calling quadprog' or not, and
> how?

Do you have a good (and feasible) solution in
advance?

If you do have one, then X0 is the argument you
would provide. You cannot do anything unless you
have one.

HTH,
John D'Errico

Subject: How do I set x0 in quadprog to get better solution?

From: William

Date: 22 Mar, 2006 05:23:47

Message: 3 of 30

er, I do provide X0 when calling quadprog, and sometimes a solution is
obtained with exitflag=1, yet such a advice is still given:' To
continue solving the problem with the current solution as the starting
point, set x0 = x before calling quadprog '.

I am afraid I'm not sure about the meaning of ' set x0 = x ' in the
above advice --- does it suggest that I do the quadprog several times
and each time use the last time's solution as this time's X0?
thank you again ... ...

Subject: How do I set x0 in quadprog to get better solution?

From: John D'Errico

Date: 22 Mar, 2006 14:08:06

Message: 4 of 30

In article <1143033827.196242.199450@i39g2000cwa.googlegroups.com>, "William" <thaenk@gmail.com> wrote:

> er, I do provide X0 when calling quadprog, and sometimes a solution is
> obtained with exitflag=1, yet such a advice is still given:' To
> continue solving the problem with the current solution as the starting
> point, set x0 = x before calling quadprog '.
>
> I am afraid I'm not sure about the meaning of ' set x0 = x ' in the
> above advice --- does it suggest that I do the quadprog several times
> and each time use the last time's solution as this time's X0?
> thank you again ... ...

Yes. I believe that is the gist of it. Restart quadprog
at the final point from the last call.

There is no assurance it will improve the solution.

John


--
The best material model of a cat is another, or preferably the same, cat.
A. Rosenblueth, Philosophy of Science, 1945

Those who can't laugh at themselves leave the job to others.
Anonymous

Subject: How do I set x0 in quadprog to get better solution?

From: William

Date: 23 Mar, 2006 05:01:08

Message: 5 of 30

Thank you for your kind & useful suggestions ~~~

Subject: How do I set x0 in quadprog to get better solution?

From: kamuran turksoy

Date: 25 Aug, 2011 21:47:50

Message: 6 of 30

John D'Errico <woodchips@rochester.rr.com> wrote in message <woodchips-B0BD65.09080622032006@syrcnyrdrs-01-ge0.nyroc.rr.com>...
> In article <1143033827.196242.199450@i39g2000cwa.googlegroups.com>, "William" <thaenk@gmail.com> wrote:
>
> > er, I do provide X0 when calling quadprog, and sometimes a solution is
> > obtained with exitflag=1, yet such a advice is still given:' To
> > continue solving the problem with the current solution as the starting
> > point, set x0 = x before calling quadprog '.
> >
> > I am afraid I'm not sure about the meaning of ' set x0 = x ' in the
> > above advice --- does it suggest that I do the quadprog several times
> > and each time use the last time's solution as this time's X0?
> > thank you again ... ...
>
> Yes. I believe that is the gist of it. Restart quadprog
> at the final point from the last call.
>
> There is no assurance it will improve the solution.
>
> John
>
>
> --
> The best material model of a cat is another, or preferably the same, cat.
> A. Rosenblueth, Philosophy of Science, 1945
>
> Those who can't laugh at themselves leave the job to others.
> Anonymous


I have a problem is kindly same this one,

In my case, i use,
[Q,val,ex]=quadprog(H_r,f_r,A_r,B_r,[],[],[],[],Q_old);

if ex==0
    options = optimset('MaxIter',5000,'LargeScale','off');
[Q,val,ex]=quadprog(H_r,f_r,A_r,B_r,[],[],[],[],Q_old,options);
end

But as you see in my codes although i use starting point, sometimes matlab gives warning as:
Maximum number of iterations exceeded; increase options.MaxIter.
To continue solving the problem with the current solution as the
starting point, set x0 = x before calling quadprog.

The thing if i already set x0=x why still it gives this error?

Regads

Subject: How do I set x0 in quadprog to get better solution?

From: Matt J

Date: 26 Aug, 2011 14:18:13

Message: 7 of 30

"kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j36fu6$k8t$1@newscl01ah.mathworks.com>...
>
> The thing if i already set x0=x why still it gives this error?
=====================

But you have not set x0=x as far as I can see. The warning (not an error) message is implying that you should do something like this

 if ex==0
    Q_old=Q;
   [Q,val,ex]=quadprog(H_r,f_r,A_r,B_r,[],[],[],[],Q_old,options);
 end

Subject: How do I set x0 in quadprog to get better solution?

From: Matt J

Date: 26 Aug, 2011 14:24:28

Message: 8 of 30

"William" <thaenk@gmail.com> wrote in message <1143032161.187275.68010@i40g2000cwc.googlegroups.com>...
> Hi everyone! I've recently used 'quadprog' function to solve an
> optimization of a 300X300 Hessian matrix:
> X=quadprog(H,f,[],[],Aeq,beq,LB,UB,X0), and Matlab 7 tends to switch to
> a medium-scale method because of the restriction conditions. Yet the
> process is slow and get a bad solution with exceeding the maximum
> number of iterations. And Matlab 7 gives advices:
> 'Maximum number of iterations exceeded; increase options.MaxIter.To
> continue solving the problem with the current solution as the starting
> point, set x0 = x before calling quadprog.'
===============

It sounds like you might have a poorly conditioned problem. Do you know what
cond(H) is, or is H too large to compute it?
 

Subject: How do I set x0 in quadprog to get better solution?

From: kamuran turksoy

Date: 26 Aug, 2011 17:20:32

Message: 9 of 30

"Matt J" wrote in message <j38aas$85t$1@newscl01ah.mathworks.com>...
> "William" <thaenk@gmail.com> wrote in message <1143032161.187275.68010@i40g2000cwc.googlegroups.com>...
> > Hi everyone! I've recently used 'quadprog' function to solve an
> > optimization of a 300X300 Hessian matrix:
> > X=quadprog(H,f,[],[],Aeq,beq,LB,UB,X0), and Matlab 7 tends to switch to
> > a medium-scale method because of the restriction conditions. Yet the
> > process is slow and get a bad solution with exceeding the maximum
> > number of iterations. And Matlab 7 gives advices:
> > 'Maximum number of iterations exceeded; increase options.MaxIter.To
> > continue solving the problem with the current solution as the starting
> > point, set x0 = x before calling quadprog.'
> ===============
>
> It sounds like you might have a poorly conditioned problem. Do you know what
> cond(H) is, or is H too large to compute it?
>

H=
  1.0e+005 *

    0.0002 0.0002 0.0000 0.0000 0.0001 -0.0000 -0.0001 0.0000 -0.0277
    0.0002 0.0003 0.0000 0.0000 0.0002 -0.0000 -0.0001 0.0000 -0.0347
    0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000
    0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000
    0.0001 0.0002 0.0000 0.0000 0.0001 -0.0000 -0.0001 0.0000 -0.0233
   -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0057
   -0.0001 -0.0001 -0.0000 -0.0000 -0.0001 0.0000 0.0001 -0.0000 0.0162
    0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0001
   -0.0277 -0.0347 -0.0000 -0.0000 -0.0233 0.0057 0.0162 -0.0001 4.8096

Here is my H in one sampling and cond(H) is 4.5950e+024.

Subject: How do I set x0 in quadprog to get better solution?

From: Matt J

Date: 26 Aug, 2011 17:52:13

Message: 10 of 30

"kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j38kl0$d5q$1@newscl01ah.mathworks.com>...
>
>
> Here is my H in one sampling and cond(H) is 4.5950e+024.
=====================

Need I say more?

Subject: How do I set x0 in quadprog to get better solution?

From: Bruno Luong

Date: 26 Aug, 2011 18:21:28

Message: 11 of 30

"Matt J" wrote in message <j38mgd$jft$1@newscl01ah.mathworks.com>...
>
>
> Need I say more?

It's not so simple to draw a rapid conclusion Matt. It is quite not straight forward the condition number of H that is important in constrained optimization. For the quadratic programming, I believe the sensitivity of the solution to errors depends the condition number of:

[H A^t;
 A 0 ]

Where A := [Aeq; Aneq(lambda==0,:)], and lambda is the Lagrange's multiplier at the solution (A is the active constrained).

Another evident that can make us to think is this: Linear programming is one particular case of Quadratic programming with cond(Hessian) = Inf ! Fortunately it still can be solved in many cases in finite-accuracy of floating points.

Bruno

Subject: How do I set x0 in quadprog to get better solution?

From: kamuran turksoy

Date: 26 Aug, 2011 19:33:10

Message: 12 of 30

"Matt J" wrote in message <j38mgd$jft$1@newscl01ah.mathworks.com>...
> "kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j38kl0$d5q$1@newscl01ah.mathworks.com>...
> >
> >
> > Here is my H in one sampling and cond(H) is 4.5950e+024.
> =====================
>
> Need I say more?

Yes you may need to say more, what is the problem? how can i over come?

Subject: How do I set x0 in quadprog to get better solution?

From: Matt J

Date: 26 Aug, 2011 20:30:33

Message: 13 of 30

"kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j38sdm$a72$1@newscl01ah.mathworks.com>...
> "Matt J" wrote in message <j38mgd$jft$1@newscl01ah.mathworks.com>...
> > "kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j38kl0$d5q$1@newscl01ah.mathworks.com>...
> > >
> > >
> > > Here is my H in one sampling and cond(H) is 4.5950e+024.
> > =====================
> >
> > Need I say more?
>
> Yes you may need to say more, what is the problem? how can i over come?
==================

First, do as Bruno suggested and apply COND to the matrix

[H A^t;
 A 0 ]

where A := [Aeq; Aneq(lambda==0,:)],
 
If it is very large, like cond(H), then your problem is ill-conditioned. Either your objective function and constraints do not contain enough modelling information to well-identify a solution. Or the quadratic function is a lot more sensitive to some variables than others, in which case you should think about rescaling some variables, expressing them in different units.

Subject: How do I set x0 in quadprog to get better solution?

From: Matt J

Date: 26 Aug, 2011 20:38:27

Message: 14 of 30

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <j38o78$pi1$1@newscl01ah.mathworks.com>...
>
> It's not so simple to draw a rapid conclusion Matt. It is quite not straight forward the condition number of H that is important in constrained optimization. For the quadratic programming, I believe the sensitivity of the solution to errors depends the condition number of:
>
> [H A^t;
> A 0 ]
>
> Where A := [Aeq; Aneq(lambda==0,:)], and lambda is the Lagrange's multiplier at the solution (A is the active constrained).
================

True, but the reported symptoms (slow and incomplete convergence) together with the huge cond(H) are very suspicious.
 
Also, if the algorithm produces a sequence where over a stretch of iterations none of the constraints are active (interior point, active set???) won't speed be influenced by cond(H) over those iterations?




> Another evident that can make us to think is this: Linear programming is one particular case of Quadratic programming with cond(Hessian) = Inf ! Fortunately it still can be solved in many cases in finite-accuracy of floating points.
========

But does the quadratic programming algorithm become similar to the linear programming simplex algorithm as cond(Hessian)--->inf?

Subject: How do I set x0 in quadprog to get better solution?

From: Bruno Luong

Date: 26 Aug, 2011 21:04:29

Message: 15 of 30

"Matt J" wrote in message <j39083$mi7$1@newscl01ah.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <j38o78$pi1$1@newscl01ah.mathworks.com>...
> >
> > It's not so simple to draw a rapid conclusion Matt. It is quite not straight forward the condition number of H that is important in constrained optimization. For the quadratic programming, I believe the sensitivity of the solution to errors depends the condition number of:
> >
> > [H A^t;
> > A 0 ]
> >
> > Where A := [Aeq; Aneq(lambda==0,:)], and lambda is the Lagrange's multiplier at the solution (A is the active constrained).
> ================
>
> True, but the reported symptoms (slow and incomplete convergence) together with the huge cond(H) are very suspicious.
>
> Also, if the algorithm produces a sequence where over a stretch of iterations none of the constraints are active (interior point, active set???) won't speed be influenced by cond(H) over those iterations?

If none of the constraints are active, then yes. Is it the case here?

>
> > Another evident that can make us to think is this: Linear programming is one particular case of Quadratic programming with cond(Hessian) = Inf ! Fortunately it still can be solved in many cases in finite-accuracy of floating points.
> ========
>
> But does the quadratic programming algorithm become similar to the linear programming simplex algorithm as cond(Hessian)--->inf?

A well designed QP algorithm must treat seriously the eventual semi-positiveness of the augmented "Hessian" [H A'; A 0], in short it must somehow deal very similar to linear programming (interior point method).

Bruno

Subject: How do I set x0 in quadprog to get better solution?

From: Matt J

Date: 26 Aug, 2011 21:43:20

Message: 16 of 30

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <j391ot$r9e$1@newscl01ah.mathworks.com>...
>
> > Also, if the algorithm produces a sequence where over a stretch of iterations none of the constraints are active (interior point, active set???) won't speed be influenced by cond(H) over those iterations?
>
> If none of the constraints are active, then yes. Is it the case here?
==============

That was part of my question. As I recall, none of the algorithms necessarily enforce constraints at every iteration, except maybe the interior point algorithm and then only the bound constraints. Even in that case, you don't know if the bound constraints are active.

Subject: How do I set x0 in quadprog to get better solution?

From: Bruno Luong

Date: 26 Aug, 2011 21:55:29

Message: 17 of 30

"Matt J" wrote in message <j3941n$4jn$1@newscl01ah.mathworks.com>...
>
>
> That was part of my question. As I recall, none of the algorithms necessarily enforce constraints at every iteration, except maybe the interior point algorithm and then only the bound constraints. Even in that case, you don't know if the bound constraints are active.

IFAIK, all algorithms keeps two or three set of active constraints, binding constraints at all iteration. The sets are not necessary updated at all iteration, but their existence always are, and they certainly updated at the final iteration.

Bruno

Subject: How do I set x0 in quadprog to get better solution?

From: kamuran turksoy

Date: 26 Aug, 2011 22:09:27

Message: 18 of 30

"Matt J" wrote in message <j39083$mi7$1@newscl01ah.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <j38o78$pi1$1@newscl01ah.mathworks.com>...
> >
> > It's not so simple to draw a rapid conclusion Matt. It is quite not straight forward the condition number of H that is important in constrained optimization. For the quadratic programming, I believe the sensitivity of the solution to errors depends the condition number of:
> >
> > [H A^t;
> > A 0 ]
> >
> > Where A := [Aeq; Aneq(lambda==0,:)], and lambda is the Lagrange's multiplier at the solution (A is the active constrained).
> ================
>
> True, but the reported symptoms (slow and incomplete convergence) together with the huge cond(H) are very suspicious.
>
> Also, if the algorithm produces a sequence where over a stretch of iterations none of the constraints are active (interior point, active set???) won't speed be influenced by cond(H) over those iterations?
>
>
>
>
> > Another evident that can make us to think is this: Linear programming is one particular case of Quadratic programming with cond(Hessian) = Inf ! Fortunately it still can be solved in many cases in finite-accuracy of floating points.
> ========
>
> But does the quadratic programming algorithm become similar to the linear programming simplex algorithm as cond(Hessian)--->inf?


As you suggest i changed units of parameters, and i am sure about objective function (there is no problem with it), i also tried unconstrained case but still when i check hessian matrix sometimes it is not positive semi definite. and cond(H) is still very large and one more strange thing to see whether hessian is positive semi definite i check eig((H+H')/2) and i see that all eigs are very small and close to each other comparing to last eig value. I could not figure out what cause this issue.

Subject: How do I set x0 in quadprog to get better solution?

From: Matt J

Date: 26 Aug, 2011 22:29:12

Message: 19 of 30

"kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j395in$8qf$1@newscl01ah.mathworks.com>...
>
> As you suggest i changed units of parameters, and i am sure about objective function (there is no problem with it), i also tried unconstrained case but still when i check hessian matrix sometimes it is not positive semi definite. and cond(H) is still very large and one more strange thing to see whether hessian is positive semi definite i check eig((H+H')/2) and i see that all eigs are very small and close to each other comparing to last eig value. I could not figure out what cause this issue.
=======================

You might want to describe the problem for us a bit more, e.g. what H means physically and how it is constructed. Same with the constraints.

Subject: How do I set x0 in quadprog to get better solution?

From: kamuran turksoy

Date: 26 Aug, 2011 23:04:10

Message: 20 of 30

"Matt J" wrote in message <j396no$bvd$1@newscl01ah.mathworks.com>...
> "kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j395in$8qf$1@newscl01ah.mathworks.com>...
> >
> > As you suggest i changed units of parameters, and i am sure about objective function (there is no problem with it), i also tried unconstrained case but still when i check hessian matrix sometimes it is not positive semi definite. and cond(H) is still very large and one more strange thing to see whether hessian is positive semi definite i check eig((H+H')/2) and i see that all eigs are very small and close to each other comparing to last eig value. I could not figure out what cause this issue.
> =======================
>
> You might want to describe the problem for us a bit more, e.g. what H means physically and how it is constructed. Same with the constraints.

I am working on system identification and using recursive identification method. Actually there are some function in matlab which do recursive least square (RLS) identification but in these function there is no constraints in parameter estimation. And in modeling and control systems parameters that i am trying to find should be between some interval. That's why i used constrained recursive identification.

Identification

My model is:

y(k)+a1*y(k-1)+a2*y(k-2)+...+an*y(k-n)=b1*u(k-1)+...+bm*u(k-m)+e(k)+c1*e(k-1)+...+cmn*e(k-mn)

here y(k) is my plant output and y(k-1) the output belongs to one step before,
u(k) input of plant and e(k) is disturbance.

And i can show my model as linear regression expression which is:

ym(k)=phi'*theta+e(k) where, ym (model output)

phi=[-y(k-1) -y(k-2) ... -y(k-n) u(k-1) ... u(k-m) e(k-1)...e(k-mn)]'
theta=[a1 a2 ... an b1 b2 ... bm c1 c2 ...cmn]'

to make my model best i have find optimal theta which makes the error between real output and my model minimum

err=y-ym(k)=y-phi'*theta

For unconstrained case RLS gives theta as (k sampling time)

theta(k)=theta(k-1)+K(k)*[y(k)-phi(k)'*theta(k-1)];
K=P(k-1)*phi(k)*[1+phi(k)'*P(k-1)*phi(t)]^-1;
P=[I-K(k)phi(k)']*P(k-1);

For case to define constraints
 My objective function is

V(theta)=[theta(k)-theta(k-1)]'*P(k-1)^-1*[theta(k)-theta(k-1)]+(y(k)-phi(k)'*theta(t))^2


and from this objective function i get

H=2*[P(k-1)^-1+phi(k)*phi(k)']
f=-2[(P(k)^-1)'*theta(k-1)+y(k)*phi(k)]

min (theta(k)'*H(k)*theta(k)+f(k)'*theta(k)

and my constraints are:

c_const=[0 0 0 0 0 0 0 0 1;0 0 0 0 0 0 0 0 -1];c_lim=[1-0.01;1-0.01];
a_const=[0 1 0 0 0 0 0 0 0;-1 -1 0 0 0 0 0 0 0;1 -1 0 0 0 0 0 0 0];a_lim=[1-0.01;1-0.01;1-0.01;];
 b_const=[0 0 1 0 0 0 0 0 0;0 0 0 1 0 0 0 0 0;0 0 0 0 1 0 0 0 0;0 0 0 0 0 1 0 0 0;0 0 0 0 0 0 1 0 0;0 0 0 0 0 0 0 1 0;0 0 -1 0 0 0 0 0 0;0 0 0 -1 0 0 0 0 0;0 0 0 0 -1 0 0 0 0;0 0 0 0 0 -1 0 0 0;0 0 0 0 0 0 -1 0 0;0 0 0 0 0 0 0 -1 0];
b_lim=[-0.00001;-0.00001;-0.00001;-0.00001;-0.00001;-0.00001;1-0.01;1-0.01;1-0.01;1-0.01;1-0.01;1-0.01];

A_r=[a_const;b_const;c_const];
B_r=[a_lim;b_lim;c_lim];

And i use quadprog as:

[Q,val,ex]=quadprog(H_r,f_r,A_r,B_r,[],[],[],[],Q_old);


I hope it is not confusing but here it a summary of codes that i use.

Subject: How do I set x0 in quadprog to get better solution?

From: Bruno Luong

Date: 27 Aug, 2011 06:23:26

Message: 21 of 30

"kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j398pa$hfd$1@newscl01ah.mathworks.com>...

>
> I hope it is not confusing but here it a summary of codes that i use.

I don't understand few details on how you derive the matrix for QP from unconstrained case. But if your derivation is correct you should get similar control vector theta (between QP and the iterative LS) if all the bounds of the inequality constraints in quadprog are relaxed to +infinity .

Have you try that to see if you coding is bug-free?

Bruno

Subject: How do I set x0 in quadprog to get better solution?

From: kamuran turksoy

Date: 27 Aug, 2011 06:42:10

Message: 22 of 30

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <j3a2gu$rc8$1@newscl01ah.mathworks.com>...
> "kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j398pa$hfd$1@newscl01ah.mathworks.com>...
>
> >
> > I hope it is not confusing but here it a summary of codes that i use.
>
> I don't understand few details on how you derive the matrix for QP from unconstrained case. But if your derivation is correct you should get similar control vector theta (between QP and the iterative LS) if all the bounds of the inequality constraints in quadprog are relaxed to +infinity .
>
> Have you try that to see if you coding is bug-free?
>
> Bruno

I guess you did not understand how to derive constraint matrix? if so,
In identification theory to have a stable model, the roots of
1+a1*q^-1+a2*q^-2+...+an*q^-n (roots respect to q) should be inside of unit circle.
In my system since i dont have negative input signal and since increase in iptup signal should decrease my output, i had to make max value for b1 b2... bm as -0.00001 and during the iteration sometimes b values go to negative infinity that's why i also add a lower constraint as -0.99 (actually i am not sure about this value it can be higher or lower) and for c parameters the same thing as a parameters.

Subject: How do I set x0 in quadprog to get better solution?

From: Bruno Luong

Date: 27 Aug, 2011 06:50:16

Message: 23 of 30


>
> I guess you did not understand how to derive constraint matrix?

No, I don't understand how you get H. But that part of work if yours, I don't even want to know how you do it.

Bruno

Subject: How do I set x0 in quadprog to get better solution?

From: kamuran turksoy

Date: 27 Aug, 2011 07:03:10

Message: 24 of 30

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <j3a438$29k$1@newscl01ah.mathworks.com>...
>
> >
> > I guess you did not understand how to derive constraint matrix?
>
> No, I don't understand how you get H. But that part of work if yours, I don't even want to know how you do it.
>
> Bruno

It is not that complicated just matrix and vector multiplication for
V(theta)=[theta(k)-theta(k-1)]'*P(k-1)^-1*[theta(k)-theta(k-1)]+(y(k)-phi(k)'*theta(t))^2
 gives you H and f

Subject: How do I set x0 in quadprog to get better solution?

From: Bruno Luong

Date: 27 Aug, 2011 07:14:10

Message: 25 of 30

"kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j3a4re$4dp$1@newscl01ah.mathworks.com>...

>
> It is not that complicated just matrix and vector multiplication for
> V(theta)=[theta(k)-theta(k-1)]'*P(k-1)^-1*[theta(k)-theta(k-1)]+(y(k)-phi(k)'*theta(t))^2
> gives you H and f

Again, I'm *not* interested in how you compute it. But have you make any check consistency between: (1) quadprog without constraint and, (2) this algo RLS:

theta(k)=theta(k-1)+K(k)*[y(k)-phi(k)'*theta(k-1)];
K=P(k-1)*phi(k)*[1+phi(k)'*P(k-1)*phi(t)]^-1;
P=[I-K(k)phi(k)']*P(k-1);

?

Bruno

Subject: How do I set x0 in quadprog to get better solution?

From: kamuran turksoy

Date: 27 Aug, 2011 07:23:09

Message: 26 of 30

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <j3a5g2$637$1@newscl01ah.mathworks.com>...
> "kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j3a4re$4dp$1@newscl01ah.mathworks.com>...
>
> >
> > It is not that complicated just matrix and vector multiplication for
> > V(theta)=[theta(k)-theta(k-1)]'*P(k-1)^-1*[theta(k)-theta(k-1)]+(y(k)-phi(k)'*theta(t))^2
> > gives you H and f
>
> Again, I'm *not* interested in how you compute it. But have you make any check consistency between: (1) quadprog without constraint and, (2) this algo RLS:
>
> theta(k)=theta(k-1)+K(k)*[y(k)-phi(k)'*theta(k-1)];
> K=P(k-1)*phi(k)*[1+phi(k)'*P(k-1)*phi(t)]^-1;
> P=[I-K(k)phi(k)']*P(k-1);
>
> ?
>
> Bruno

Yes, i checked and they gave exactly same results.

Subject: How do I set x0 in quadprog to get better solution?

From: Bruno Luong

Date: 27 Aug, 2011 07:35:27

Message: 27 of 30

"kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j3a60t$7hj$1@newscl01ah.mathworks.com>...

>
> Yes, i checked and they gave exactly same results.

That means quadprog converge without constraint?

Then I don't understand why the constraints introduced degrade the convergence of quadprog. It should improve IMO. Unless the recursive formula falls down [ again there is a question how to make right the correspondence between RLS and QP recursion ] and you run into a non-obvious issue that in turn causes the hessian to become singular during the recursion.

Bruno

Subject: How do I set x0 in quadprog to get better solution?

From: Matt J

Date: 27 Aug, 2011 15:33:10

Message: 28 of 30

"kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j398pa$hfd$1@newscl01ah.mathworks.com>...
>
> My model is:
>
> y(k)+a1*y(k-1)+a2*y(k-2)+...+an*y(k-n)=b1*u(k-1)+...+bm*u(k-m)+e(k)+c1*e(k-1)+...+cmn*e(k-mn)
>
>[SNIP]
>
> to make my model best i have find optimal theta which makes the error between real output and my model minimum
>
> err=y-ym(k)=y-phi'*theta
==============

This looks wrong. The model equation as you first wrote it is equivalent to
phi'*theta=0. So it looks like err should really be

err=phi'*theta

In any case, how many time samples do you use?



> and my constraints are:
>
> c_const=[0 0 0 0 0 0 0 0 1;0 0 0 0 0 0 0 0 -1];c_lim=[1-0.01;1-0.01];
> a_const=[0 1 0 0 0 0 0 0 0;-1 -1 0 0 0 0 0 0 0;1 -1 0 0 0 0 0 0 0];a_lim=[1-0.01;1-0.01;1-0.01;];
> b_const=[0 0 1 0 0 0 0 0 0;0 0 0 1 0 0 0 0 0;0 0 0 0 1 0 0 0 0;0 0 0 0 0 1 0 0 0;0 0 0 0 0 0 1 0 0;0 0 0 0 0 0 0 1 0;0 0 -1 0 0 0 0 0 0;0 0 0 -1 0 0 0 0 0;0 0 0 0 -1 0 0 0 0;0 0 0 0 0 -1 0 0 0;0 0 0 0 0 0 -1 0 0;0 0 0 0 0 0 0 -1 0];
> b_lim=[-0.00001;-0.00001;-0.00001;-0.00001;-0.00001;-0.00001;1-0.01;1-0.01;1-0.01;1-0.01;1-0.01;1-0.01];
>
> A_r=[a_const;b_const;c_const];
> B_r=[a_lim;b_lim;c_lim];
==================================

Note - most of these constraints are actually bound constraints. For faster performance, you should really be using quadprog's LB and UB input arguments to specify them.

Subject: How do I set x0 in quadprog to get better solution?

From: kamuran turksoy

Date: 27 Aug, 2011 19:37:11

Message: 29 of 30

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <j3a6nv$98a$1@newscl01ah.mathworks.com>...
> "kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j3a60t$7hj$1@newscl01ah.mathworks.com>...
>
> >
> > Yes, i checked and they gave exactly same results.
>
> That means quadprog converge without constraint?
>
> Then I don't understand why the constraints introduced degrade the convergence of quadprog. It should improve IMO. Unless the recursive formula falls down [ again there is a question how to make right the correspondence between RLS and QP recursion ] and you run into a non-obvious issue that in turn causes the hessian to become singular during the recursion.
>
> Bruno

Yes, it converge, but the problem is, without constraint it can give parameters which makes model unstable which can not be used in control and identification

Subject: How do I set x0 in quadprog to get better solution?

From: Matt J

Date: 28 Aug, 2011 22:46:09

Message: 30 of 30

"kamuran turksoy" <kamuranturksoy@gmail.com> wrote in message <j398pa$hfd$1@newscl01ah.mathworks.com>...
>
> to make my model best i have find optimal theta which makes the error between real output and my model minimum
>
> err=y-ym(k)=y-phi'*theta
>
> For unconstrained case RLS gives theta as (k sampling time)
>
> theta(k)=theta(k-1)+K(k)*[y(k)-phi(k)'*theta(k-1)];
> K=P(k-1)*phi(k)*[1+phi(k)'*P(k-1)*phi(t)]^-1;
> P=[I-K(k)phi(k)']*P(k-1);
>
> For case to define constraints
> My objective function is
>
> V(theta)=[theta(k)-theta(k-1)]'*P(k-1)^-1*[theta(k)-theta(k-1)]+(y(k)-phi(k)'*theta(t))^2
===================

It looks to me as though P(k) and K(k) come from Kalman filtering formulas, but it's not entirely clear to me what motivates adding this term

[theta(k)-theta(k-1)]'*P(k-1)^-1*[theta(k)-theta(k-1)]

to your error term
y(k)-phi(k)'*theta(t))^2

Is it for regularization? If so, shouldn't you be applying some weighting to the regularization term to control its relative contribution to V(theta)?

Furthermore, do you even know whether P(k-1) is stably invertible?
What is cond(P(k-1))?

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us